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Preface 



This course aims to make students aware of the physical origins of the main partial differential 
equations of classical mathematical physics, including the fundamental equations of fluid and 
solid mechanics, thermodynamics, and classical electrodynamics. These equations form the 
backbone of modern engineering and many of the sciences, and solving them numerically is 
a central topic in scientific computation. 

The goal of these notes is to establish the essential ideas with maximum clarity and 
reasonable rigor, with a conventional approach to notation and nomenclature, to serve as 
an introduction to the enormous literature on these topics. I have not attempted a compre- 
hensive treatment of, say, continuum mechanics - such a treatment at the upper division - 
beginning graduate level of this course would occupy more than the semester available for 
this topic in our curriculum. Many excellent comprehensive texts in continuum mechanics, 
as well as in classical thermodynamics and in electromagnetism, will be accessible to the 
student who finishes this course. My goal is rather to get across the core idea cutting across 
all of these subjects - that the fundamental field equations result from combining a few ba- 
sic physical principles, some clearly stated assummptions and approximations, and vector 
calculus. 

I have also not presented this subject as a branch of modern analysis, with the rigor that 
would entail - although such an exposition is certainly possible [Marsden and Hughes, 1994] . 
Moreover, the information I have developed here about the solutions of these equations is 
quite limited. The big three examples in PDE - the Laplace, wave, and heat equations - turn 
up as special cases, and I do establish a few basic facts about their solutions. Additional 
properties of solutions, and methods to approximate them numerically, are developed in 
other parts of the CAAM and MATH curricula: 

• CAAM 437, MATH 423: analysis of PDEs; 

• CAAM 452, CAAM 552: numerical methods for PDEs and their analysis. 

Students should enter this course with at least a nodding acquaintance with vector cal- 
culus, as taught in the fourth semester of a calculus sequence, and some background in dif- 
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ferential equations, equivalent to that obtained from CAAM 336. For applied mathematics 
students to gain optimally from the study of this material, mastery of rigorous multivariate 
advanced calculus (as provided by succcessful completion of a course similar to CAAM 402) 
is essential. 

Nothing in my exposition of this fundamental, well-aged material is novel. I have bene- 
fitted greatly from the excellent text of M. Gurtin on continuum mechanics [Gurtin, 1981], 
and drawn liberally from many other references, citations to which are sprinkled throughout 
the text. In general, I have resisted the temptation to "improve" on the notational con- 
ventions common in the subject. My principal aim is to prepare students in this course for 
further reading and participation in the very active fields of science and engineering that 
base themselves on this material. While nothing like a standard syntax for topics of these 
notes exists, I have chosen notation that is as close in spirit to the "norm" as I could come. 

William W. Symes 
Houston, Texas 
August 2006 

Added January 2012: this version of the notes differs from the previous in that I have made 
a few minor changes of notation, added some problems, and enlarged the trailing chapters 
in various ways. The aim and spirit of the enterprise remains precisely the same. 
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Chapter 1 



Kinematics and Conservation 



Several of the physical principles underpinning continuum mechanics take the form of con- 
servation laws, which state that some property of the material in motion is preserved during 
the motion, or changes in a prescribed way. The fundamental postulate of continuum me- 
chanics - that the properties of matter can be localized in arbitrarily small volumes - has 
its mathematical expression in the use of densities, i.e. amount per unit volume, to repre- 
sent mass, momentum, energy, etc. Materials compress, expand, and change shape during a 
motion, so the change in a property density also has a purely kinematic component, simply 
due to the motion. This chapter is devoted to the expression of basic conservation principles 
as partial differential equations, which account for both physical and geometric causes of 
density changes. 



1.1 Deformation and Motion 



This section introduces some of the mathematics necessary to describe motion of 3D con- 
tinua. These are mostly standard ideas from vector calculus, expressed in perhaps unfamiliar 
notation suited to the application. 

A "material body" is the volume occupied by a specific piece of material. Note that this 
definition says both more and less than it appears to say. Namely, it implies some notion of 
what a "material" is, and how it may be divided into pieces - without, however, supplying a 
clear description of either idea. In fact, the real content of the definition is simply that pieces 
of material, whatever that is, are identified with the subsets of space which they occupy. A 
further, fundamental notion is that any given piece of material may occupy (at different 
times, or under different conditions) many different subsets of space. These various subsets 
are the possible configurations of the body. One of the great ideas of this subject is that one 
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of these possible configurations, that is, subsets of space, should be singled out as a reference 
configuration. 

In this discussion, "space" is Euclidean 3-space, or R 3 in common notation. Say Q C 
R 3 - is the reference configuration of a material body. We use capital letters to denote 
coordinates in this reference configuration or location. In effect, these coordinate vectors 
X = (Xi,X2,X3) T label positions within the body. Because some material is supposed to 
be located at each position, these vectors XgQ are called material points. So, once again, 
a (vague) physical notion, the material located at a point in space, is confounded with a 
precise geometrical one, namely the point's coordinate vector. 

Another brilliant idea, which came about sometime in the development of this subject, is 
that the various possible configurations of a material body must be related by mappings, or 
vector-valued functions: if uj is a configuration of the body whose reference configuration is 
Q, then each point x£w must be the position occupied by the material which is located at 
some other position X G Q in the reference configuration. That is, if it's the same material, 
each bit ("point") of it must be in both configurations, and it must be possible to say, 
unambiguously, that x 6 w is the same material ("point") as X G fl. That is, there must 
be a mapping x : — > R 3 so that x(f2) = u, and any material point X G fi represents 
the same material as is located at x(X) G uj. Such a map is called a deformation, as the 
"new" configuration uj of the material is presumably a deformed version of the reference 
configuration Q. 

I draw your attention to several peculiarities of the notation. I've used capital letters, 
as noted before, to single out vectors X and their coordinates X±,X 2 ,X 3 in the reference 
configuration, and small letters x = (xi, x 2 , x%) T for vectors describing points in arbitrary 
(deformed) configurations. This notational device is standard in the subject - see for ex- 
ample [Marsden and Hughes, 1994, Chung, 2007]. Also conventional, at least implicitly, is 
the confusion of vector variables and vector functions: depending on context, x can denote 
either the vector variable x = (xi, x 2 , x^) T or the vector-valued function x : Q — > R 3 . This 
unfortunate habit at least has the virtue of reducing the amount of notation: Marsden and 
Hughes, for example, write instead <fi : Q — > R 3 for the deformation function, but implicitly 
identify with x whenever convenient. By simply using x as the name of the deformation 
function as well as a typical coordinate vector in the deformed configuration, I avoid in- 
troducing additional notation like Marsden-Hughes' 0, which soon disappears anyway, but 
confront you with the necessity of determining meaning from context. This will be good 
practice for reading the literature of continuum mechanics. Finally, I insist that vectors 
are columns, unless otherwise indicated - that is, x = (xi, x 2 , xs) T , not = (xi,X2,xs). This 
insistence is not conventional, but has an excellent justification, to appear shortly. 

Deformations must satisfy several more conditions, which express natural physical no- 
tions in geometrically precise ways. Deformations must be invertible: physically, it must be 
possible to reverse a deformation to recover the reference configuration from a deformed one. 
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Figure 1.1: A deformation map x and its inverse X. 

Once again confusing variables and function names, I will use X : uj — > R 3 as the name of the 
inverse of x : Q — > R 3 , that is, X = x _1 . Thus if Y G f2, then X(x(Y)) = Y. Deformations 
(and their inverses) must be continuous: this is a precise statement of the notion that if 
points X, X' are close, then the corresponding points x(X),x(X') in the deformed material 
should also be close - no ripping or tearing allowed! See Figure 1.1 

Even more is required: until further notice, deformations are assumed (often tacitly) to 
be continuously differentiable, of class C 1 , with inverses also of class C 1 . 

The derivative of x, identified with its Jacobian matrix of partial derivatives, is conven- 
tionally denoted by Vx in this subject: 

Ox ■ 

(Vx)ij = i,j = 1,2,3 

That is, the element is the infinitesimal rate of change of the ith component of deformed 
position vector with respect to the jth component of the position vector of the same bit of 
material ( "material point" ) in the undeformed configuration. 

Note that many other notations (Dx, J x ,...) are used in other contexts for the same 
thing. However, the V notation is universal in continuum mechanics. 

The derivative of a deformation function x is called the deformation gradient - note that 
this object already has a perfectly good name ("derivative"), but gets a special one for use 
just in continuum mechanics! 

I'll use the same notation for arbitrary vector-valued functions of vector arguments, with 
arbitrary numbers of components in domain and range: that is, V denotes the derivative, 
expressed as a matrix relative to some coordinate system, and the ith row contains the 
partial derivatives of the ith component, whereas the jth column contains the jth partial 
derivatives of all of the components. 

With this notation, the chain rule takes a simple form: if u : R m — > R p (p-vector function 
of an m- vector), and v : R™ — > R m (m- vector function of an n- vector) are differentiable, 
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then the composition u(v(x)) is also differentiable, and 

V(u(v(x))) = (Vu)(v(x))Vv(x) 

i.e. the derivative of the composition of two functions is just the matrix product of the 
derivatives of the two components in the order of composition. 

Note that u and v don't need to be defined in all of the respective Euclidean spaces for 
the composition be well defined and for the chain rule to hold. It's enough that the range of 
v be contained in the domain of u, so that the composition makes sense. 

[In particular, if X : (a, b) — > Q is differentiable curve defined on the interval (a, b) C R, 
then 

|(x(X(i))=Vx(X(i))^(i). 

This works because coordinate vectors like X are columns - that's why it's the right conven- 
tion!] 

Notational Remark: A more satisfactory way to write this relationship suppresses the 
role of the dummy argument (x above). The composition of u and v is denoted u o v: that 
is, for any x in the domain of v, u o v(x) = u(v(x)). In this more natural notation, which 
gives a name to the composition, the chain rule reads 

V(uov) = (Vuov)Vv 

Unfortunately this more rational approach to notation for continuum mechanics is confined 
to a few authors - for example [Marsden and Hughes, 1994]. The vast bulk of the literature 
on continuum mechanics uses the older standard syntax, with dummy variables and indices 
galore. Since the purpose of these notes is to prepare you, dear reader, for the literature 
as it actually is, rather than as the author wishes it were, I will stick with the notational 
conventions introduced in this chapter - less than ideal, but representative of what you will 
find in papers, books, and presentations on this subject. 

There is one annoying exception to this nice organization of partial derivatives into arrays. 
For a scalar function $(x), the symbol V$(x) denotes the gradient, which is the column 
vector of partial derivatives. That is, the gradient is an object of the same type as a position 
vector. This is not the natural specialization of the vector V notation: a scalar function 
would naturally be identified with a 1-vector-valued function, and its derivative represented 
as a row vector. A row vector is the matrix representation of a linear map, from R 3 to R. 
The derivative of a real-valued function is exactly such a linear map, so indeed the natural 
representation of the derivative should in principle be a row vector. However the transpose 
of this row vector is the Riesz representer of the derivative: its inner product with a direction 
(vector) w e R 3 gives the directional derivative in that direction, 

uT $(x + ew)-f(x) 
(V$(x)) J w = lim^o— — 
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For various reasons, the Riesz representer (the fastest increase direction of in the Eu- 
clidean metric) is more useful than the natural matrix representation of the derivative, and 
is conventionally signified by V (or grad). 

So V$ is the transpose of what you would expect, and the chain rule reads 

V($(v(x)) = Vv(x) T (V$)(v(x)). 

You should satisfy yourself that this matrix product actually does express the compo- 
nentwise chain rule 

^ (v(x ))^» (v(x) )g (x ) 

I'll also need the change of variables formula for multiple integrals: if u is the image of 
Q under the deformation x, then 

/ cfe$(x)= f dX$(x(X))detVx(X) 

[Note that on the LHS of this equation, x is a dummy variable representing an arbitrary 
point in u, whereas on the RHS it refers to the deformation (vector- valued function) mapping 
Q to u, and X is the dummy variable. This type of contextual reading will occur over and 
over again.] 

Finally, I will often use the following Averaging Theorem: if / and g are (scalar or vector 
valued) continuous functions defined on the same ("reasonable") domain f2 C R 3 , and for 
any cube C ClVl 

I f=f 9 

Jc Jc 

then /(X) = g(X) for all XgO. 

Note that if you divide both sides by the volume of ui, the hypothesis becomes: the 
averages of / and g over any cube are the same - thus the name of the theorem. 

Other classes of sets, other than cubes, for example all balls, or even more general classes 
of sets, will work as well. 

The meaning of "reasonable" here, i.e. the sets over which integration has a well-defined 
meaning, is a complicated issue which I am not going to address in this course. Cubes, 
balls, images of these under continuously differentiable and invertible maps, and unions and 
intersections of finite numbers of these, are certainly "reasonable". Courses on measure 
and integration, such as MATH 425 at Rice, make these notions precise. Good texts are 
[Rudin, 1986, Royden, 1988, LANG, 1993, ?]. 
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One final condition on deformations: they are required to be orientation preserving, 
which means that 

detVx > 

This condition prevents deformation from turning the material "inside out" . An example of 
a non-orientation preserving map figures in one of the exercises - according to this condition, 
it cannot define a deformation. 

A lot of this subject will concern motions of material continua. A motion of a material 
body is simply a time-parameterized family of deformations, x(-,t), defined for an interval 
(to,ti) of "time" (R!) - that is, for each t G (t ,ti). x(-,i) : Q — > R 3 is a deformation. The 
range, or image, of the reference configuration under the motion is a time-parameterized 
family uj t = x(f2,t) of subsets of R 3 , that is, a curve in the set of configurations. 

Motions will be assumed different iable (continuously, if necessary several times) without 
further mention, as functions of all arguments. 



1.2 Conservation of Mass, Material Description 

For present purposes, mass is simply an attribute of material bodies (that is, subsets of 
R 3 ) with certain properties, to be enunciated. Most of the properties flow from assuming 
the existence of a mass density function p rcf : Q — > R on the reference configuration. This 
function is assumed to be continuously differentiable and positive. It gives the mass of any 
"reasonable" part pCfl through 

m (p) = / dX p Tci (X). 

It follows from what we have assumed so far that there also exists a mass density for the 
deformed states, p(x, t). The integral of which over a part of the deformed body gives the 
mass of that part: for p C ou t , 

/ dx p(x, t) 

is the mass in p at time t. 

This is particularly true if p = p t , where p t = x(p, t) = {x(X, t) : X G p} is a moving 
part of the body, i.e. the volume occupied at time t by the material which occupied some 
part p of the undeformed body. See Figure 1.2. 

VERY IMPORTANT: Throughout this section p t will denote such a moving part - i.e. 
it is a time-parameterized family of subsets of R 3 corresponding to a part of the material, 
moving with the material. The same material is in p t for all t. In particular p t is not a fixed 
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Figure 1.2: A body in the undeformed state (Q) and a part of it (p), and at two times t < ti 
during a motion. The moving part p t , or t — to,ti and any other t, is presumed to contain 
exactly the same material as the corresponding part p in the undeformed configuration of 
the body. 

subset of R 3 except in case the motion preserves it, and its variation with t is determined 
by the motion and is not arbitrary. 

So we can state a (material) description of the conservation of mass as follows: 

/ dX p re f(X) = / dxp(x,t) 

Jp Jpt 

The change-of-variable rule for multiple integrals gives an alternate form in terms on inte- 
gration over p only: 

/ dXp Ici (X)= [ dXdetVx(X,t)p(x(X,t),t) 
Jp Jp 

This relation must hold for "any" subset pcfl. Certainly it holds when p is a sufficiently 
small cube. The Averaging Theorem then implies that the (continuous) integrands must 
actually be identical: 

p rcf (X) = detVx(X,f)p(x(X,t),t). 
Since the x(-, t) is invertible and its Jacobian is everywhere positive, this implies in turn that 

p(x,t) = p rcf (X(x,t))detVX(x,t). 

[Recall that X is the inverse deformation, i.e. the inverse function to x.] 

The preceding equation is the material description of the conservation-of-mass principle: 
it defines a quantity (mass density) in the deformed configuration explicitly in terms of the 
reference configuration and the relation between the two. It is also possible to produce a 
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spatial description of this conservation law, that is, one that refers only to quantities defined 
in the deformed configuration (in this case, p(x, t), not mentioning the reference configuration 
at all. The principle is general, and is explained in the next section. 



1.3 Reynolds' Transport Theorem 



Reynolds' Transport Theorem provides a general rule for passing between material and 
spatial descriptions of conservation principles. 

Suppose that x : x (t ,ti) — > R 3 is a motion, p t = x(jo, t) C oo t = x(f2,t) is a moving 
part, and <3>(x, t) is a continuously differentiable function on an open set in R 4 containing 
{x, t : x G p t , t < t < ti}. Then the function 



Q(t) = I dx®(x,t) 
Jpt 



is differentiable. The main result of this section is an expression for its rate of change in 
terms of partial derivatives of $ and of the motion. 

The easiest way to proceed is to refer the integration back to the undeformed configura- 
tion, so that at least the domain does not depend on t: 



Q(t)= [ rfXdetVx(X,t)$(x(X,t),t) 

J V 



The t-derivative will then simply pass under the integral sign: 



dQ 
dt 



(t) = J p dX ^detVx(X,t)) $(x(X,t),t) + detVx(X,t) (JU(x(X, 



Evidently we will need to understand how to differentiate determinants of matrix-valued 
functions of t. In fact 

^-detVxfX, f) = detVx(X,t)tr f ^-Vx(X, t)(Vx(X, t))' 1 ) , 
at \at J 

where tr denotes the trace, i.e. the sum of the diagonal entries. I defer the proof of this 
formula to the next section. 

Since x and X are inverse functions, their Jacobians are inverse matrices as noted above, 

so 

= detVx(X,t)tr fv x |(X,t)(VX)(x(X^),^ 



The chain rule says that for any vector-valued function w, 

Vw(X(x,f)) = (Vw)(X(x,t))VX(x,t) 

so substituting x = x(X, t) on both sides, 

(Vw(X(x,t))) x=x( x, t ) = Vw(X)(VX)(x(X,t),t) 

The right hand side of the last expression for the time derivative of the Jacobian determinant 
has the form of the right hand side above, with w = dx/dt. Thus 



JUetVx(X,t) = det Vx(X, t)trV f ^(X(x,t),t) 



x=x(X,t) 



An interesting piece of this expression has a significance of its own: define the spatial 
velocity field: 

v(x,*) = ^(X(x,i),f). 

This is actually the velocity of the "material point" which happens to be located at x at time 
t. Indeed, the velocity of the point corresponding to location X in the reference configuration 
is 

V(X,t) = ^(X,t); 

this field is known as the material velocity. If a "material point" is at x at time t, then it 
was at X(x, t) in the reference configuration, so the velocity is precisely v(x, t) given above. 
That is, 

v(x,t) = V(X(x,t),f). 

Recall from vector calculus the divergence operator V-: for any 3- vector field w, 

3 d Wl „ 

v • w = L = trVw - 

i=l ux i 

With these notations, the expression for the time derivative of the Jacobian determinant 
turns into the remarkable formula 

d 

— detVx(X, t) = detVx(X, i)(V • v)(x(X, t),t) 

This handles the first half of the expression for the time derivative of Q. The second half 
involves 

r) fhc 8<t> 

-$(x(X,t),t) = — (X,t) T (V$)(x(X,t),0 + ^(x(X,t),t) 

= (^ + v-V$) (x(X,t),t) 



This combination of partial derivatives occurs so often that it gets a special name and 
notation: it is called the material derivative, represents the derivative along the trajectory 
x(X, t) of a material particle, and is denoted 

DO „ 
— = — + v • V 

Dt dt 

Putting all of this together, 

^§(t) = J rfXdetVx(X,t) + (V • v)*) (x(X,f),f) 

This is a setup for changing the integration back to p t . This leads to the final expression for 
the rate of change of Q, which is Reynolds' Transport Theorem: 

d f r { D<& \ 



L ix (f 



+ V-(v») (x,i) 



The last equality follows from one of the vector calculus identities which you will establish 
in Problem Set 1. 



1.4 Differentiating Determinants 

Suppose that A(t) is a continuously differentiable n x n matrix valued function of t, defined 
on an open interval containing t = 0. In this section I'll derive a surprisingly simple formula 
for the derivative of det A(t) at t — 0. This formula played an important role in the derivation 
of Reynolds' Transport Theorem. 

Before establishing determinant derivative formula in general, I should mention a special 
case in which you can see the result without any trouble at all. Suppose that A(t) is 
continuously differentiable and diagonalizable, with continuously differentiable diagonalizing 
matrix T(t) so that 

T-\t)A(t)T(t) = diag(ai(t),...,a n (t)) 
with continuously differentiable and positive eigenvalues a±, ...,a n . Then 

det A(t) = ai(t)a 2 (t)...a n (t) 

so 

^-detA(t) = ^(t)a 2 (t)...a n (t) 
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+... + a 1 (t)a 2 (t)...^(t) 

. , . / 1 da.\ , . 1 da n , \ 

= detA{t) . . , (t) + ... + .. , t 

= deU(f)tr ^(t)A-\t)y 

In fact the formula just derived is true, regardless of whether A is diagonalizable or not 
(though without that assumption, it requires a different proof!). 

It is necessary to assume that detA(O) > 0, from which it follows that det A(t) is positive 
for small t. Note that for our application, in which A(t) = Vx(X, t), this is one of the 
defining properties of a motion. 

The first order Taylor expansion of A is 

A(t) = A(0)+tA'(0) + R(t) 

where R(t) = o(t). Thus 

deU(t) = deU(O) det(J + t(A- 1 (0)A'(0) + fU-^O)^))) 

Set A = —1/t. The scaling property of determinants (det(cA) = c n detA) gives 

detA(t) = (-t) n detA(0)det(A/ - A- 1 (0)(A'(0) + r 1 ^*))) 

Fact from Linear Algebra: the coefficient of A n_1 in the characteristic polynomial det (XI — 
M) of an n x n matrix M is the negative of the trace trM = Ma. So 

detA(t) = (-t) n detA(0)(A n - A n " 1 trA(0)- 1 (A'(0) + t" 1 ^*)) + P(A) 

where P(X) represents a polynomial in A of degree < n — 2. Multiplying through bt (— t) n = 
A"". 

= deti4(0)(l + ttrA(0)- 1 (A , (0) + t^Rtf)) + 0(f)) 
= deU(0)(l + ttrA(0) _1 j4'(0) + o(f)) 
Comparing with the Taylor series 

. . . c/det^4 . . . . 
= detA(O) + *— — (0) + o(t) 

(JjL 

and noting that there is nothing special about t — gives the desired identity: 

— - — = detAtrM" 1 — ) 
rft y dt 1 
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(since tiAB = tiBA for any n x n matrices A, B). Otherwise put, 

^-logdetA = trM" 1 ^) 
dt v dt 1 



1.5 Conservation of Mass, Spatial Description 

This follows directly from the Transport Theorem. Simply take $(x, t) = p(x, t) - from the 
material description of conservation of mass, 



Q(t) = P= Pref 

Jp t Jp 



is independent of time, being the mass "transported" in the deforming part p t of the material. 
So 



This must hold for "all" (regular enough?) parts p t - it certainly holds when p t is a small 
cube. That is enough (Averaging Theorem!) to show that 

= l + v ' (v " ) = ir + (v '^ 

which gives two forms of mass conservation in the spatial description. 



1.6 Volume change and the significance of the velocity 
divergence 

You get the volume of a region in R 3 by integrating the constant function 1 over it. So for 
our moving part p t , 

V(t) = I dxl 

You can apply Reynolds' Transport Theorem to this integral too, as $ = 1 is a perfectly 
legitimate choice of function. Of course there is not reason for the volume to remain constant, 
and indeed moving parts may expand or contract during a motion. The rate of change of 
volume is 

^-(t) = [ dxV ■ v(x,t). 
dt Jp t 
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That is, the velocity divergence is the density for volume rate of change! 

With this insight, you can see the meaning of the terms in the RHS of the Transport 
Theorem: the first (D<fr/Dtdx) is the intrinsic rate of change of the density (&dx) itself 
along material particle trajectories, whereas the second ((V • v)&dx) reflects the rate of 
density change due solely to changing geometry, i.e. contraction or expansion of the volume 
occupied by the material. 

A very important class of motions also reveals its characterisation courtesy of this identity. 
A motion is volume preserving if and only if for any part p t , the volume is constant in time: 
dV/dt = 0. The identity above and the Averaging Theorem imply that this is the case 
exactly when 

V • v(x,t) = 

A material is called incompressible if it admits only volume-preserving motions. Many com- 
mon fluids, eg. water, are approximately incompressible under ordinary "room" conditions. 



1.7 Conservation of Linear Momentum 



We don't need to assume the existence of a momentum density - this is assured by the 
existence of a mass density, and even uniquely determined. The mass per unit volume at x 
at time t is p(x, t), and the velocity of material at that time and place is v(x, t). Therefore 
the momentum density, i.e. momentum per unit volume, of the material at x at time t is 
p(x, i)v(x, t), and so the momentum of a part p t of the body at time t is 

PpM = / dx p(x,t)v(x,t) 
Jpt 



Note that this is already a spatial description, not referring explicitly to the location of 
material particles in the undeformed body. As an excercise, find a material description of 
momentum. 

Momentum is not an invariant of motion - instead, its evolution is governed by Newton's 
second law of motion, the famous "F=ma", interpreted as applying to a material body or 
part of it. "Mass times acceleration" is the same as "rate of change of momentum". Thus 
Newton's law for the part pt of the material reads 



dt 



(t) = F\p t 



where the right hand side is simply shorthand for the sum of all forces acting on the part p t 
at time t. These forces will get careful attention shortly. Meanwhile express this relation in 
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terms of p and v, using the vector version of Reynolds' Transport Theorem (Problem 1-9): 

D 

dx 



F[ Pt ] = J dxj^(pv) + (V-v)pv 



This statement uses the material derivative applied to a vector field (the momentum density 
pv: for any vector field (vector-valued function) w, 

Dw dw .„ . 
_ s _ + (Vw)v 

An alternate statement of the vector Transport Theorem is also possible, as in the scalar 
case. 

If the force acting on p t were given by integration of a density, f say, over p t , then you 
could use the Averaging Theorem: since p t is arbitrary, you would conclude that 

f = ^(pv) + (V-v)pv 

In Problem Set 1 you will show that a consequence of mass conservation is the identity 

_ (p4 >) + (V.v)p* = / ,— 

which holds for any scalar- or vector- valued function $. Thus momentum balance takes the 
quite pleasing form 

f Dv 

f = P Dt 

which is about as close to "F=ma" as one ought to expect to get. 

This is the form that most of the famous PDEs ("laws of motion") take, such as the 
Euler and Navier-Stokes equations and the various versions of elastodynamics. It will take 
some work to understand the internal structure of the force density f . 

In fact it's actually possible to view momentum conservation in a way that looks even 
more like classical mechanics. The center of mass of the body is the average position within 
the body, weighted (as it were!) by the mass density: 



q(i) = — [ rfxxp(x, t) 
m Jut 



>u>t 

where 

m — dx p(x, t) 

JU} t 

is the mass of the object. Note that we have integrated over the entire body, occupying 
volume oo t at time t, so no mass is lost or gained - and the m in the last equation is really 
as independent of t as the notation indicates!. 
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The rate of change of the center of mass is 



dq. 1 r £)x 

at m Ju) t Ut 



According to the definition of the material derivative for vector-valued functions (like x 
here a dummy variable denoting position, not the material particle trajectory x(X, t)\), 

Dx <9x .„ . 

_ = _ + (Vx)v = v 

since the Jacobian of the position coordinate vector x is I. Thus 

■(t) = — / dx p(x, i)v(x, t) 

Tfl Jut 



~dt 



Recalling the definition of momentum p(t) (in this case p t = u> t ) and clearing the denomi- 
nator, you get 

dq 

m Tt =v 

which reads "mass times velocity equals momentum". Thus from this point of view, the 
body behaves as if it were a classical mechanical particle with mass m concentrated at the 
center of mass, whose momentum is the total momentum of the body. In this notation, 
conservation of momentum amounts to 

dp d 2 q 
F — — — m-— 

dt dt 2 

which really is "F=ma" ! 
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1.8 Problems 



Problem 1-1: Suppose that fl (the volume occupied by the undeformed body) is the unit 
cubic meter: {X : Om < x 0ji < lm, i — 1,2, 3}, which deforms by stretching uniformly under 
the motion 

x(X,t) = (l + t)X. 

(a uniform expansion or dilatation). The initial density is homogeneous, p re f = 1000kg/m 3 . 
(a) Verify that this formula defines a motion for t > 0. Compute (b) v(x, t) (units of 
m/s), (c) p(x,t). 



Problem 1-2: Define / : R 3 ->■ R 3 by 

/(x) 

and set 



£o = 



a;? - 2rr 3 

3 1 

2s| 



2 
1 

-1 



Find the best linear approximation (i.e. the first order Taylor series) of / at X. 



Problem 1-3: Define the map x : B 3 x R — > R 3 by 

x(X,t) = 



Xq : \ + tXy! 
X0.2 + t(x _ 2 + 2X , 3 ) 
^0,3 + ^0,1 



where B 3 = {X e R 3 : ||X|| < 1} is the unit ball in R 3 . 

(i) Show that x is a motion for small t. Is x a motion for all t? For all positive t? [Hint: 
The issue is whether x is invertible on the reference configuration, i.e. with its given domain. 
If X e B 3 , then | ^o,i | < 1- If t is small enough, the "other" root of the quadratic giving Xo,i 
in terms of x\ is too big to be the first component of a vector in the unit ball.] 



(ii) Define <p : R 3 x R ->■ R by 



0(x, t) = 2x\ + x 2 + tx 3 



Compute (a) (<90/di)(x(X, t), t), (b) d/dt((j)(x.(X,t),t)) (i.e. formulas, good for all X, £ - 
note carefully the order of the parentheses! In one case you compute the partial derivative, 
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then substitute x = x(X, £), whereas in the other you substitute first then compute the 
partial derivative). 

(iii) What would be involved in computing v(x, t) for this motion? Can you give a formula 
for v? ["No" is an acceptable answer, if you can justify it - otherwise produce v(x, t)\] 



Problem 1-4: Take Q to be all of R 3 except the origin, and define 

X 



x X 



X 



Note that this map sends points inside the unit sphere to points outside, and points near 
the origin to points very far away - it turns Q inside out! Show that 

detVx(X) = -||X||- 6 

and conclude that this map is not a deformation, even though its Jacobian determinant does 
not vanish anywhere in its domain - it is not orientation preserving! 



Problem 1-5: Suppose that the motion x(X, t) is volume-preserving, and that the density 
at time t — is independent of position x (homogeneous), i.e. p(x, 0) = constant (indepen- 
dent of x). Show that the density remains homogenous at all times, i.e. p(x, t) = constant 
(independent of x and t). Hint: Use the spatial description of the conservation of mass to 
conclude that Dp/Dt = 0. Then interpret the material derivative as rate of change along 
the trajectory of a material particle. 



Problem 1-6: Suppose that $(x, t) is a differentiable scalar- or vector-valued function, 
defined on a subset of R 4 including all positions of the body. Show that for any part p t , 

— J dxp{x,t)®{x,t) = dx p(x,i) — (x,t) 

Deduce that 

D D$ 
-(p*) + (V.v)p* = p— . 



Problem 1-7: Suppose that Q (the volume occupied by the undeformed body) is the 
cylinder: {X : 0m < x 0t3 < lm, Xq j + x\ 2 < 1}, and that the body undergoes the rotation 
about the xo,3 axis with angular frequency u: 

x\ (X, t) = xo,i cos out + xq.2 sin out 
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x 2 (X, t) = —x 0j i sin out + x 0j2 cos out 

X 3 (X,t) = X ,3 

The initial density is homogeneous, p ref = 1000kg/m 3 . (a) Verify that this formula defines 
a motion for all t. Compute (b) v(x, t) (units of m/s), (c) V ■ v(x, t) (what units?), (d) 
p(x,i). 



Problem 1-8: Derive the following vector calculus identities, from the definition of gra- 
dient and divergence. All, when reduced to components, are straightforward consequences 
of the Leibnitz rule for partial derivatives of scalar- valued functions. In each statement, / 
and g denote scalar-valued functions and u and v (3-)vector-valued functions. Recall that 
the gradient is defined (inconsistently!) for a scalar- valued function / as a column vector 
with ith row entry (V/)i = J^, and for a vector-valued function u as the 3x3 matrix 

m x dui 

(Vu) - = W- 



The divergence of a vector-valued function u is 

3 

— -\ 

dx 



3 dv 

i=i 



I'll use both u ■ v and u T v to denote the standard Euclidean inner product J2 u i v i of vectors 
u,v. 

Here are the identities you need to justify: 

1. V(fg) = fVg + gVf 

2. V-(/u) = V/-u + /V-u 

3. V(u • v) = (Vu) T v + (Vv) T u 

4. V(/u) = u(V/) T + /Vu 



Problem 1-9: Derive the vector version of Reynold's Transport Theorem: for a vector- 
valued function w(x, t), and a moving part p t within its domain of definition, 

d f , f , Dw 

dt jp t jp t 

in which the material derivative Dw/Dt is defined by 



f f Dw 

/ dxw = dx — — h (V • v)w 

Jpt Jpt Dt 



Dw dw .„ , 
s + (V.), 
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Chapter 2 



Strain 



Deformations which don't change the internal distances in a body are in some sense trivial: 
they should not give rise to changes in the internal forces which each part of the body exerts 
on its complement. Thus there is need for a construction which expresses not deformation 
per se, but the extent to which a deformation or motion gives rise to changes in the internal 
geometry of a body. 



2.1 Rigid Deformation and the Definition of Strain 

The official adjective for a deformation which does not change internal distances is rigid. To 
say that x : il — > R 3 does not change internal distances means 

||x(X 1 )-x(X 2 )|| = HXi-Xall 

where the vertical double bars denote Euclidean distance. A rigid motion is of course a 
motion consisting of rigid deformations. 

Suppose that the two points under consideration are close, say X and X + AX. The 
Taylor series of x gives 

x(X + AX) - x(X) = Vx(X) AX + R 

with R = o( || AX ||). So the change in length (squared) is 

||x(X + AX) - x(X)|| 2 - ||AX|| 2 

= AX T (Vx(X) T Vx(X) - J)AX + o(|| AXf) 
= 2AX T L(X)AX + o(||AXf) 
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where L(K) is the (Lagrangian or) material strain tensor 

L(X) = i(Vx(X) T Vx(X)-7) 

As with everything else in this subject you can ask where the stuff currently at x came from, 
instead of where the stuff originally at X went to. That is, express the change of length as 

||Ax|| 2 - ||X(x + Ax) -X(x)|| 2 = 2Ax T e(x)Ax + o(||Ax|| 2 ) 

where the (Eulerian or) spatial strain tensor is 

£(x) = i(/-VX(x) T VX(x)) 
The following result expresses the relation between the strain tensors and rigidity: 

Theorem Suppose that the interior Q° is non-empty. Then the following conditions are 
equivalent: 

1 . x : — > R 3 is a rigid deformation 
2. 

x(X) = x(Y) + M (X - Y) 

where Y G f2° and M is a rotation matrix, i.e. an orthogonal matrix with positive 
determinant, M T M = MM T = I, detM = 1 - that is, x is a Euclidean affine trans- 
formation; 

3. The Lagrangian strain tensor L vanishes; 

4. The Eulerian strain tensor E vanishes. 



Proof: 1 =>- 2: Pick an arbitrary Y G Then the hypothesis means that for any XgO°, 

(x(Y) - x(X)f (x(Y) - x(X)) - (Y - X) T (Y - X). 

Both sides are differentiate in a neighborhood of Y; a calculation gives the derivatives (with 
respect to Y): 

Vx(Y) T (x(Y) - x(X)) = Y - X 

Since x is a deformation, Vx(Y) T is invertible. Write M = (x(Y) T ) _1 . Then solving the 
above equation for x(X) gives 

x(X) = x(Y) + M(X — Y). 
20 



Invoking the rigidity condition again, for any Xsfl, 



(M (X - Y)) T (M(X - Y)) = (X - Y) T (M T M)(X - Y) 
= (X- Y) T (X- Y) 

Since Y is interior to Q, this equation holds for X in a sufficiently small ball centered at Y. 
Since any vector is a multiple of X — Y for some X in such a ball, it follows that M is an 
orthogonal matrix: M T M = I. Since / is a deformation, its gradient must have positive 
determinant, therefore det M > 0. 

2 ^> 3: the deformation gradient Vx(X) is constant and = M, so, 

L = M T M -1 = 

because M is orthogonal. You can also easily see that VX(x) = M T , so E = 0. q.e.d. 

Thus the strain tensors really do measure the amount by which a deformation deviates 
from rigid motion - this deviation is strain. Other quantities expressing changes in the 
internal state of the body, such as internal forces, should depend on the strain, not on the 
deformation. 

2.2 Strain and Displacement 

It's also convenient to express strain in terms of the displacement of points during the de- 
formation, defined as 

w(X) = x(X) - X 

(Lagrangian displacement) or 

u(x) = x — X(x) 
(Eulerian displacement). The expressions are 

L = i(Vw + Vw T + Vw T Vw) 

E = -(Vu + Vu T - Vu T Vu) 

Homogeneous deformations are described by 

x(X) = x a + A(X - x 6 ) 

for some choice of vectors x a ,Xb and a matrix A e R 3 . Since Vx(X) = A, this formula 
defines a deformation if det A > 0. Then 

X(x) = X a + A- 1 (x-X 6 ) 
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(What are X a ,X b in terms of x a ,xt?) It's easy to compute the strains: 

L = \(A T A-I), E= 1 -(I-(A-rA-') 

where I've left off X and x because both of these matrix-valued functions are constant! 

The Polar Decomposition Theorem of linear algebra states that for any square matrix A, 
there are symmetric matrices S and T and an orthogonal matrix R of the same size as A so 
that 

A — RS — TR. 

If detA > then R is a rotation matrix, and detS, .detT > also. 

Substitution in the formulae for the strains and use of the identities satisfied by rotation 
matrices results in the formulae 

L = \{S 2 -I\ E= 1 -(I-T- 2 ) 

Similarly, for an arbitrary deformation you can write 

Vx(X) = R(X)S(X) = T(X)R(X) 

in which the symmetric and rotation parts of the polar decomposition are now matrix-valued 
functions. The same manipulations lead to 

L(X) = i(S(X) 2 - /), E(x) = \{I- T(q(x))- 2 ) 

That is: the strain tensors depend only on the symmetric parts of the polar decomposition 
of the deformation gradient Vx(X) . 

Since the strains are symmetric matrices, the Spectral Theorem of linear algebra asserts 
the existence of an orthonormal basis of eigenvectors, called the principal axes of Eulerian or 
Lagrangian strain. The continuum mechanics name for the eigenvalues is principal strains. 
Some of the problems in the next set should give you a sense of the meanings of principal 
strains and axes. 



2.3 Problems 
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Problem 2-1. A uniform extension in the direction e e R 3 is a homogeneous deformation 
of the form 

x = x a + A(X - x 6 ), A = 7 + (A - l)ee T 

where A is a positive scalar. Suppose that e = (1, 0, 0) T . Compute the strain tensors L and 
E explicitly (i.e. give the matrix entries in terms of the scalar A), and compute the principal 
strains and axes (eigenvalues, eigenvectors of L, E). Draw the unit cube {X : < x 0ji < 1} 
and its image under a uniform extension (do your best a perspective drawing) - can you see 
why this deformation is called an extension? 



Problem 2-2. A homogeneous deformation with 



A 



1 7 
1 
1 



where 7 is a scalar, is called pure shear in the (1,2) (or (x,y)) plane. Compute L, E, and 
the principal strains and axes. Draw the image of the unit cube under this deformation. Is 
"shear" the right word? 



Problem 2-2 - Extra Credit. Conjecture a relation between Eulerian and Lagrangian 
principal strains and axes, based on the result of the last two problems. Prove it. 



Problem 2-3. Suppose that / : R -> R, f(0) = 1, and 

x(X,t) = 



1 + f 
f(t) 
1 



X 



Specify (i.e. give a formula for) the function f(t) so that this motion is volume-preserving. 
Describe the motion thus defined in words, as best you can. 



Problem 2-4. Suppose that A is any 3x3 matrix, i.e. A E R 3x3 . Explain why the 
formula 

x(X,t) = x a + (/ + L4)(X-x 6 ) 
defines a motion for small t. Compute i>(x, t). 



Problem 2-4 - Extra Credit. Under what conditions on the matrix A is the motion 
defined above volume-preserving? 
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Problem 2-5. Suppose that : f2 — > R 3 is any continuously differentiable function on the 
undeformed material volume, and define 

$(x,i) = 0(q(x,t))detVq(x,t) 

Show that 

— (x,t) + (V-t;)(x,t)$(x,t) =0 
(N.B.: this is a review problem, so it really is as easy as it looks!) 
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Chapter 3 



Stress 



3.1 Zoology of forces 

The forces acting on continuum materials divide logically into those exerted in the material 
by itself, and those which are the result of some external agency. For a part occupying 
subvolume u of the material in Q, we can write 



As was explained in Ch. 1, momentum balance takes the form of a differential equation if 
the force can be expressed as an integral of a force density over u, so it is interesting to see 
when that is the case. 

I will have little to say about external forces: I will simply assume that they are given in 
terms of a force density f having units of force per unit volume: 



Newtonian gravity is an important example of such a force; near the surface of the Earth 
the gravitational force density is —pge r , in which e r is the "radial" direction, i.e. the unit 
vector (in the units used to define p and g) on the line through the point on the surface and 
the Earth's center, pointing outward. 

The standard analysis of internal forces, as presented in texts like [Gurtin, 1981], depends 
on two fundamental hypotheses, both asserting the locality of internal forces, i.e. the absence 
(or unimportance!) of action at a distance: 

Internal Forces are Contact Forces: There is a contact force density function t[oj] 
on the boundary du, having units of force per unit area, giving the force exerted by the 



FtotalM = Fintfw] + F cxt [uj] 
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remainder of the body (material in Q 



uj) on the material in ou: 




I tL4(x)t[w](x) 



Here cL4(x) is the area element on du. 

This equation expresses the conviction that all influence of one part of the material on 
the remainder takes place at the its boundary - what is going on in the remainder away from 
the boundary has absolutely no influence in the force experienced by the part. 



Cauchy's Hypothesis on the Form of Contact Forces: The contact force density 
depends only on the tangent plane, or equivalently the (outward) normal vector n(x) at 
each point of the boundary dui: 



for a function T : VL x R 3 — > R 3 , called the traction. 

This apparently innocent hypothesis has profound consequences. Note that whereas the 
contact force density of the first hypothesis depends on the part uj in some apparently very 
general way, Cauchy's hypothesis ties it to a much more restricted form: the contact force at 
a point on the boundary x e du does not depend on the shape or curvature of the boundary 
at x but only on its orientation, ie. on the direction of the normal at x. This implies that 
all parts whose boundaries are tangent at x experience the same contact force density at 
x. This is clearly a very great restriction in the form of contact forces - just how great will 
become apparent shortly. 

The intuitive meaning of Cauchy's hypothesis is a little harder to divine, but it also ex- 
presses locality of influence. For purposes of computing contact force, it is as if at each point 
of the boundary, one replaces the boundary by a small flat piece of the tangent plane. If 
the curvature or higher order shape descriptors were involved, then the contact force would 
depend less locally on the disposition of material. Since the orientation of the boundary is 
clearly a factor in the contact force experienced by real materials (think of sticking your 
hand out the window of a speeding car - the force exerted by the airstream depends dra- 
matically on its orientation!), Cauchy's hypothesis asserts the maximum locality of contact 
forces consistent with experience. 



Remark: In all of this discussion, it has gone without saying that the boundary du is 
sufficiently regular to have a well-defined normal. I will make this assumption throughout. 



tM(x)=T(x,n(x)) 
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3.2 Covariance of forces 



Up to now, the discussion has treated the coordinate system as given. Of course coordinates 
are more or less arbitrary, and "physics" should not depend on the choice of coordinates. 
These notes will make use of only a part of the invariance of classical mechanics: all of the 
laws stated so far, and those to come, must be unchanged in their meaning under Euclidean 
changes of coordinates, i.e. new coordinates x' related to old coordinates x by 

x' = x'(x) = x' + i?(x - x ) 

where R is an orthogonal matrix. 

This is indeed a coordinate transformation, i.e. it is invertible with smooth inverse: the 
inverse Euclidean transformation is 

x = x(x') = x + i? T (x - x(,) 

Note the similarity in form, but not in meaning, to the definition of rigid deformation: 
here no deformation occurs, i.e. "the same stuff is in the same place" - just the coordinate 
description of physical location has changed. 

This difference in meaning has an important consequence: it is not necessary to require 
that R is a rotation matrix, i.e. has a positive determinant. It is perfectly acceptable to use 
R with deti? = —1, as occurs for instance changing from a right-handed to a left-handed 
coordinate system. In particular R = — / is fine (note that —I is not a rotation in 3D!). 

The velocity vectors of a material point in the two coordinate descriptions are related: 

<9x' <9x 

~dt = at 

so 

v'(x',t) = Rv(x,t) 

Note what's going on here. The velocity vector v does not simply compose with the change 
of coordinate - i.e. the velocity in the new coordinate system is not v(x', t'). Instead the 
value of this vector quantity also rotates, as described by the above equation, and in the 
same way that the position vector x — x rotates into x' — x' . A vector quantity which 
changes in this systematic way is covariant. 

A covariant vector field (vector- valued function) is not just a vector field - it is one vector 
field for each coordinate system, i.e. a vector-valued function on the Cartesian product of 
the Euclidean transformations and R 3 . Moreover, all of these vector- valued functions are 
related, as just described for the velocity. This is a subtle idea, and the notation does not 
support it well. 
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The Jacobian matrix of x' with respect to x is simply R, so the chain rule for a vector- 
valued function <fr reads 

V x *(x'(x)) = (V**)(x'(x))i2 
If $ is a covariant vector field, then 

$'(x'(x)) = i2*(x), *(x) = R T &(x!(x)) 



so 



whence 



V x $(x) = V x (i? J *'(x'(x))) = (it> J V x ,<&'R>)(x'(x)) 
at 

= i? T — + -R V x '<fr -RR f = i? T 

That is: the material derivative of a covariant vector field is also a covariant vector field. 

This applies in particular to the velocity field. So the left-hand side of the momentum 
balance law is also covariant (since p is scalar, and changes simply by composition with the 
change of coordinates): 

Dv / ,Dv'\ 



Therefore the momentum balance law can only hold regardless of coordinate system if 
force densities are covariant vector fields. 

To see how this rule applies to traction, you need to understand how normal vectors 
change under changes of coordinate. Suppose that P = {5x. : n ■ 5x = 0} describes the 
tangent plane to a surface at some position: i.e. n is the normal there. That is, each 5x e P 
is the tangent vector of some curve lying inside the surface. Then the tangent plane to the 
same surface in the new (primed) coordinates is P' = {5x.' : Sx.' = R5x forsome <5x e P}, as 
you can see by looking at the coordinates of the curves. To be perpindicular to all of these 
vectors, n' must satisfy n' = i?n, i.e. n = Rn', so normal vectors are also covariant. 

Therefore traction, being a force (density), must satisfy 

T(x, 6) = R T T'(x', 6') = R T T'(x', R6) 

i.e. 

T'(x', 9) = i?T(x, R T 6) 

for any unit vector 9. 
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Curious example: In particular, taking R = —I, x = x we get 

T(xo,0) = -T'(xo,-0) 

This looks like Newton's Law of Action and Reaction - i.e. the force exerted by the material 
outside a surface on the material inside is the same in magnitude as, and the opposite in 
sign to, the force exerted by the material inside the surface on the material outside - but 
isn't quite, since the right hand side is the traction field in a different coordinate system. 
Newton's Law in the case of tractions will follow from momentum balance, as a consequence 
of Cauchy's Theorem in the next section. 

The careful reader will have noticed that the traction field appears only with unit vectors 
(namely, outward unit normals) in its second argument, in the statement of momentum 
balance. Therefore we are free to define its values for non-unit vectors in any way we see fit. 
The following convention is convenient: 

Scaling Convention for Traction: For any real a > 0, and 9 e R 3 : 

T(x,a0) = aT(x,0) 
i.e. tractions are homogeneous of degree 1 in their second arguments. 

3.3 Cauchy's Theorem 

Cauchy's Theorem on the Linearity of Traction: There exists a 3 x 3 matrix valued 
function S'(x), called the Stress Tensor (or just stress) so that 

T(xo,0) = 5(x o )0, x en, 6 G R 3 

Proof: This is Cauchy's famous "tetrahedron argument", see [Gurtin, 1981] pp. 97-106 for 
example. It uses a number of facts about multidimensional integrals which are established 
in advanced calculus. 

The Scaling Convention, explained in the last subsection, implies that it is sufficient to 
consider only unit vectors 9. Since the traction (like everything else, whenever convenience 
dictates) is continuous at least, it also suffices to consider only 6 with nonzero components, 
the other cases following by continuity. After a Euclidean change of coordinates with diagonal 
matrix leaving x fixed: 

x = x + i?(x' - x ), R u = -f-r, R i:j = 0, 1 < i, j < 3 
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9' = R9 in the new coordinates has positive components. Since T transforms linearly, as 
estabished in the last section, the general case will follow from this one. 

Construct a tetrahedron T with one face in the plane {x + 5x. : 5x. ■ 9 = 0} (i.e. with 
normal parallel to 9) and the other three faces parallel to the three coordinate directions. 
As it is no trouble to do so, I'll assume that the deformation in question is part of a motion, 
though I will suppress t from the notation. This covers static deformations, which are 
snapshots of trivial motions (for which rates of change, including velocities, are zero). 

The balance of linear momentum states that 



where the first term on the right hand side represents the (internal) force exerted by the 
rest of the body on the stuff in the tetrahedron, and the second term represents all external 
forces. 

For each R > there is exactly one such tetrahedron for which the smallest circumscribed 
sphere centered at x has radius R; denote this tetrahedron by T R . Essentially R is the 
diameter of T R . 

Assume as always that all functions appearing in the integrals above are at least once 
continuously differentiable. Then the volume integrals are bounded by a constant times the 
volume of the tetrahedron, which is proportional to R 3 . The standard notation to express 
this fact concisely is that the volume integrals are 0(R 3 ). 

By Taylor's Theorem, 



Since an area integral of a continuous function over the bounded surface dT R is bounded by 
the area of dT R times the max of the integrand, 



The boundary dT R is the union of four faces: Fjj, perpindicular to 9, and F l R) i = 1,2,3, 
perpindicular to ei = (1, 0, 0) T , e 2 = (0, 1, 0) T , and e 3 = (0, 0, 1) T respectively. In particular 
the (outward!) normal is constant on each face: on F R , n(x) = 9, and on F R) i = 1,2,3, 
n(x) = — e,. So 





T(x, n) = T(x , n) + 0(|x - x |) = T(x , n) + O(R) 
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= A° R T(xo, 9) + A^T(xq, -d) + ^T(x , -e 2 ) + A|T(x , -e 3 ) + 0(R 3 ) 
where A R is the area of face F l R . Since A R = R 2 A\, this is the same as 

= i? 2 (A°T(x , 9) + A}T(x , -ei) + ^T(x , -e 2 ) + ^T(x , -e 3 )) + 0(R 3 ) 

Plug this expression for the integral of traction over the boundary of into the momentum 
balance law, use the already- noted fact that the volume integrals are all 0(R 3 ), and divide 
by R 2 to get 

A°T(x , 9) + A}T(x , - ei ) + ^T(x , -e 2 ) + A?T(x , -e 3 ) = 0(R) 
The left-hand side is independent of R, so letting R — > you obtain the identity 
A?T(x , 0) + A}T(x , - ei ) + A?T(x , -e 2 ) + A?T(x , -e 3 ) = 

Claim: 

-^q — 6i, i — 1, 2, 3 
The easiest proof of this fact uses the Divergence Theorem: 



Divergence Theorem For any reasonable domain Q and vector field $ defined on and 
continuously differentiable up to the boundary, 



/ diV'$(x)= / <L4(x)$(x) • n(x) 
Jn Jan 



Proof of Claim: Set <&(x) = e, in the Divergence Theorem applied to Q = r±. The 
left-hand side vanishes, while the right hand side is A\9i + A}, which proves the Claim. 

Applying this area identity to the traction identity above gives 

T(x , 9) = -0iT(xo, -ei) - # 2 T(x , -e 2 ) - # 3 T(x , -e 3 ) 



where S(x )j,i = -T(x , -e^-. 

Since the proof has made no special assumption about the coordinate system in which 
traction has been described, the above identity holds for any coordinate system and any 
vector 9 with positive coordinates in that system. 
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Consider now an arbitrary vector 9 with non-zero components in some coordinate system. 
As mentioned at the beginning of the proof, introduce new coordinates via the Euclidean 
transformation 

x' = x + R(x — x ), R = diag 

so that x'(x ) = x and 9' = R T 9 = R9 has positive components. From the covariance of 
traction, 

T(x , 9) = i? r T'(x , 9') 

= -R T (j:(9' ■ e;)T'(x , -ej) 

Note that ite; = e-, i.e. R carries the ith unit basis vector of the old coordinate system into 
the ith unit basis vector of the new coordinate system. Therefore #'-e- = R9-Rei = 9-ei = 9, h 

so 

T(x , 9) = -R(£{9 ■ e t )R T T(x , -Re 1 ,) 

i=l 

= -X)0 i T(xo,-e i ) = S , (xo)0 
i=i 

where S'(xo) is as defined above. That is, we have established the conclusion for all vectors 
9 with nonzero components. 

Finally, any vector can be approximated arbitrarily closely by vectors with only nonzero 
components. Since both sides of 

T(x ,9) = S(x )9 

are continuous in 9, this identity extends by continuity to all 9. This observation completes 
the proof, q.e.d. 

From the covariance property of traction, discussed in the last section, you can immedi- 
ately read off 

Change of Stress Tensor under change of coordinates : if coordinate systems x, x' 
are related by a Euclidean change of coordinates 

x' = Xq + i?(x - x ) 

then the stress tensors (describing the same state of internal forces) for the two coordinate 
systems are related as follows: 

S'(x') = RS(x)R T 
Another immediate consequence of Cauchy's Theorem is 
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Newton's Law of Action and Reaction for traction: 

T(xo,-0) = -T(xo,0) 

which follows immediately from the linearity of T. This observation also gives an alternate 
expression for the ith column of 5 , (x ): 

S , (xo). )i = T(x ,e i ) 

It also combines with the observation at the end of the last section to reveal a curious 
fact: if the coordinate systems x and x' are related by a Euclidean transformation with 
R = —I, i.e. the axes of the one are the axes of the other, but reversed, and x = Xq, then 

T(x o ,0) = T'(x o ,0) 

Note that this is a very special property of the coordinate-reversal transform, and does not 
hold for general Euclidean transformations. 

Remember that T(x , 6) is the contact force at a point on (any!) boundary with outward 
normal 9. Thus the traditional interpretation of the stress tensor: S(x.o)ij is the ith compo- 
nent of force exerted at x by the material in {xj > xoj} on the material in {xj < x j}, i.e. 
across the j coordinate plane. Thus the diagonal components of S are often called normal, 
the off-diagonal components tangential. 



3.4 Differential form of momentum balance 



Cauchy's Theorem also gives an alternate expression for the contact force term in the mo- 
mentum balance law, in terms of the stress tensor S: 



/ cL4(x) T(x, n(x)) = / <L4(x) S(x) n(x) 

J du J du 



I du J du 

The ith. component of this vector can be written 

/ dA(x) (S(x) v )n(x) = / dA(x) (S(x)£) • n(x) 

Jdu Jdu 

which has the same form as the right hand side of the Divergence Theorem. 

The Divergence Theorem is so important in the remainder of these notes that I will state 
it again: 
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Divergence Theorem For any reasonable domain and vector field <fr defined on f2 and 
continuously differentiable up to the boundary, 



j diV'$(x)= / dA(x)&(x) • n(x) 



So 



/ cL4(x) T(x, n(x)) i . = / dx V • (S(x) v ) T 
Now define the divergence of a matrix valued function such as S by the formula 

Then the contact force term becomes 

/ <L4(x) T(x, n(x)) = / cfeV-S(x) 

This is just a special case of the 



Divergence Theorem for Matrix Valued Functions: For any reasonable domain Q 
and matrix valued <fr defined on Q and continuously differentiable up to the boundary, 

/ & V'$(x)= / cL4(x)$(x)n(x) 
Jn Jan 

Of course the definition of divergence for matrix valued functions is set up exactly to 
make this form of the divergence theorem follow from the usual one. 

I have now expressed all three terms in the momentum balance law as volume integrals 
over the (arbitrary) part u> of the material body. The usual application of the Averaging 
Theorem gives the differential form of momentum balance: 




= V • S + f . 



This very important equation will serve as the equation of motion for most of the continuum 
theories presented in the next few chapters. 



3.5 Conservation of angular momentum and the sym- 
metry of stress. 

To define angular momentum, I require another relic from vector calculus, the cross product 
u x w of two vectors u and w: 

(u x w)i = u 2 w 3 - u 3 w 2 
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(u X w) 2 = U 3 Wi — U1W3 

(u x w) 3 = UiW 2 - u 2 w l 

Another way to describe the cross product: form the 3x2 matrix [uw]. Then the ith 
component of the cross product is the ith 2x2 minor determinant of this matrix (i.e. leave 
out the ith row, take the determinant of the remaining 2x2 matrix, and multiply by ( — 
Yet another way: form the skew-symmetric matrix 



U 



-u 3 u 2 
u 3 —ui 
—u 2 Ui 



Then u x w = Uw. 

Problem Set 3 asks you to work out a few properties of the cross product. 
The angular momentum of a moving part u> t is 

o wt (0 = / dxp(x, t)(x. x v(x, t)) 

Jut 



The external torque is 



The internal torque is 



= / dx~x. x (p(x, i)v(x, t)) 

Ju>t 

/ dxsc x f (x, t) 

Ju>t 

J (M(x)xxT(x,t,n(x,t)) 



Iduit 

where I have made the notation reflect that both T and n depend explicitly on t. 

The principle of Conservation of Angular Momentum states that the rate of change of 
angular momentum is the sum of all torques acting on the body, that is, 

^- ( dxp(x, £)x x v(x, t) — ( (L4(x)x x T(x, t, n(x, t)) + f dixxf(x,i) 

Ctt J uit J duit J u>t 

This principle is actually independent of the Conservation of (Linear) Momentum, and in 
fact the two together imply the symmetry of the stress tensor. To see how this works, first 
analyse the left hand side of the above equation: from Problem 2 of Problem Set 2, 



/ (ii/)(x,t)xxv(x,t)= [ dx p(x, t) — (x x v(x, t)) 

Ju>t JuJt L)t 

= J dxp(x,t) (xx^) 
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as follows from the vector calculus identities you will establish in Problem Set 3, especially 
Problem 4. 

On the other hand, 

f (M(x)xxT(x,t,n(x,t))= / dA(x)xxS(x,i)n(x,t) 

J duit J Quit 

= [ dxW ■ (x x S(x,t)) 

Jut 

r 3 d 
= / dx E^— ( xx s (x,t)-j) 

Ju) t j—i OXj 



= [ dx,j:C 



o ; X .S(X./).. ; • X X 



= / di(W(x)+xx(V'S(x,t)) 
where W x = S 32 - S 23 , W 2 = S 13 - S 31 , and W 3 = S 2 i - S 12 . 

Putting this all together, the law of angular momentum balance becomes 

f dxxx ( / p(x,t)^(x,t)) = f dx(W + x x (V • S(x,t) + f(x,i))) 

Juit \ J-Jt / Jut 

The Averaging Theorem, applied in the usual way, plus a little rearrangement of terms, 
shows that 

/ Dv 

W (x, t) = x x p— (x, t) - V • 5(x, i) - f (x, i) 



However the differential form of the linear momentum balance law derived in the last section 
shows that the right hand side of the above expression vanishes. In view of the definition of 
W, we have proved 



Theorem (Symmetry of Stress): S = S T . 



3.6 Problems. 



Problem 3-1: Show that for any vectors u and v, and scalar a, 
• a. u x (av) = (au) x v = a(u x v) 
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• b. u x v = —v x u 

• c. v x v = 

• d. For a matrix A e R 3x3 , denote by u x A the matrix whose jth column is u x A.j, 
i.e. the columnwise cross product. Show that 

U X v = (u x J)v 

and that the matrix U defined in the text is in fact u x /. 



Problem 3-2: Assume that u and v are continuously differentiable vector fields. Show 
that the jth column of V(u x v) is 

v(uxv) - = a^ xv + ux a^ 

= (Vu). J xv + ux(Vv). J 



Problem 3-3: With notation as in Problem 2, show that 

V(x xu) = -ux/ + xx Vu 



Problem 3-4: In this problem, v denotes the velocity field, not an arbitrary vector field. 
Use the results of Problems 1 - 3 to show that 

D Dv 

X X V = X X 

DV ' Dt 



Problem 3-5: Show that 



a. 



• b. 



x 



V(V|X|) = R J "b^ XxT 



c. (in 3D only!) 



V- = 0, x^O 
x 



V • (|x|x) = 4|x| 
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Problem 3-6: In this problem, u and v are continuously different iable vector fields, as 
in Problem 2, and Vx denotes the curl operator. 

/ \ T 

_ / duz dv>2 du\ du% dii2 du\ \ 

V x u = 



8x2 0x3 ' 8x3 dxi ' dx\ 8x2 



Show that 



• a. 

V • (u x v) = (V x u) • v — (V x v) ■ u 

• b. If u is twice continuously differentiate, then 

V x (V x u) = V(V • u) - V • Vu 
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Chapter 4 

Ideal Fluids and the Euler Equations 



4.1 Characterization of Ideal Fluids 

The characteristic property of fluids is that stress is scalar, i.e. 

S(x,t) = -p(x,t)I 
where p is a scalar function, the pressure. The minus sign is a convention. 
The expression of stress divergence in terms of p is simple: 

(v-s) z = j:^-(-pi lj ) = (-vp) t 

j=l OX j 

i.e. 

V • S = -Vp 

An ideal fluid is characterized by two further properties: 

• admits only volume preserving deformations, i.e. is incompressible: 

V- v = 

• has constant density: 

p(x,t) = po 
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These two hypotheses together imply the conservation of mass (verify this!), so the only 
independent conservation law that we have discussed so far is the conservation of momentum, 
which reads 



which reveals the nonlinearity of Euler's equations. In fact these are quite difficult to solve - 
to understand the solutions mathematically, or to construct good numerical approximations. 
Therefore a great deal of effort has been devoted to analysis of more tractable but still useful 
special cases. 



Motion of a material (fluid or otherwise) is steady if all quantities describing it are indepen- 
dent of time. For a fluid, this means that v, p, and p are all independent of time. Euler's 
equations become 

Vv v + Vp = f , V • v = 

and the first of these shows that / must be independent of time as well. 

One can learn a few things about flow without actually solving the equations of motion. 
For example, suppose that steady flow takes place in a pipe, i.e. a region Q for which the 
boundary divides into three parts, representing inlet, outlet, and sides: dQ = ri n ur out ur si d e . 
Since the sides are non-porous, the velocity vector must be tangent to the sides, and so 
v ■ n = at every point on the side r si d c . 

Since p is independent of t, the conservation of mass in spatial terms is just 



Integrating this identity over the domain of steady motion throws away lots of information, 
though not all: the divergence theorem gives 



The last three equations together are Euler's equations for an ideal fluid. 



Note that the momentum balance law can also be written 




4.2 Steady Flow 



V • (p v) = 
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(the side integral vanishing because v is perpindicular to n there). These two terms give 
the net rate of mass flow in respectively out of the pipe (actually the first term is minus the 
inflow), and conservation of mass says that they must cancel: in steady flow, inflow must 
equal outflow. 

Another interesting consequence emerges if you make another assumption that is satisfied 
approximately by flow of ordinary fluids through pipes: on Tin, v ~ — Vi n n- m and p ~ p m ; 
on r out , v ~ — v out n out an d P — Pout- Here the subscripted values of v and p are scalar 
constants, independent of position. Denote by A in and A out the areas of the inlet and outlet 
respectively. Then the mass balance law reduces to 

= -4outP0^out - AnPo^in 

SO 

^out^out = ^m^in 

which is a phenomenon anyone who has played with a garden hose will recognize: if you 
make the outlet area smaller than the inlet area, the outlet velocity will be larger than the 
inlet velocity. 

The balance of linear momentum implies that 

° = In p °( Vv ) v + Vp 
This seems intractable, but the following fact gets us out of trouble: 

Lemma: For a volume preserving motion, the velocity term in the material derivative is a 
divergence: 

(Vv)v = V • (vv T ) 

Proof: Leibnitz says 

V • (vv T ) = (Vv)v + v(V • v) 
and the second summand on the right vanishes for volume preserving motion, q.e.d. 

Thus the divergence theorem applies: 

= / p (v ■ n)v +pn 
Jan 



Note: This really is the divergence theorem, for matrix valued functions: if is a scalar 
function, 

[ V = / V • (0/) = / 0n 
Jn Jn Jen 
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Then, since the flow is tangent to the pipe walls (v • n = 0) 

0=(7 +/ +/ ) (p (v-n)v+pn) 

= A>ut(Pout + Po^out) n out + -4in(Pin ~ Po^in) n m + / 

J r s idc 

The third term has the interpretation of total force exerted by the fluid on the pipe walls. 
This expression shows that you can compute this total force merely from knowing the density, 
pressures, and velocities at the inlet and outlet. 



4.3 Circulation and Irrotational Flow 



A closed path C in a domain Q e R 3 is the image of a continuously differentiable map 
c : [0, 1] — > Q for which c(0) = c(l). The circulation of a vector field u : Q — > R 3 is the 
integral of the tangential component of u around C: 



/«.&=/ *u(c( t )).-( t ) 



IC 

The change of coordinates rule for integrals shows that this definition is independent of 
orientation preserving reparametrization of C: any two parametrizations c(t) and c'(t') 
yield the same circulation, so long as dt/dt' > 0. That is, the circulation depends only on 
the vector field, the closed curve, and the sense in which the close curve is traversed (its 
orientation). 

A motion is irrotational if the circulation of the velocity field about any closed curve in 
the domain of the motion is zero. 

A closed curve C C Q is the boundary of a smooth bounded surface D iff 



1. D is the image of a continuously differentiable injective map d from the unit disk 
U = {(xi, x 2 ) T : x\ + x\ < 1} C R 2 into Q, and 

2. c(t) = d(cos27rt, sin27rt) is a parametrization of C. 



There are two continuous unit vector fields defined on such a smooth bounded surface which 
are perpindicular at every point to its tangent field. Given the orientation of the bounding 
curve implicit in a parametrization, as in the definition, the positive normal n is the one for 
which the determinant 

, \dc d 2 c 
d6t ^'^' n 
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is positive at every point where the first two columns do not vanish. In terms of the map 
mentioned in the definition, this choice of normal is a positive multiple of 

dd dd 

To interpret irrotational motion, recall Stoke 's Theorem: if the closed curve C bounds a 
smooth piece of surface D, and n is the positive unit normal field with respect to the chosen 
orientation of C, then for any continuously differentiable vector field u defined on and near 
D, 

/ u-dx= / cL4(x)(V x u) • n 
Jc Jd 

It is an advanced calculus excercise to show how the averaging theorem implies that 
irrotational flow satisfies 

V x v = 

This is the local or differential interpretation of zero circulation. 



4.4 Potential Flow 

Statement of the next major result result requires more nomenclature. 

A homotopy h in a domain Q C R 3 is a continuous map from the unit square [0, 1] x [0, 1] 
into fl. The end curves of a homotopy are the images of the maps h , hi : [0, 1] — >■ Q defined 
by 

ho(t) = h(0,f), hi(t) = h(l,t) 

In fact, for general s G [0, 1], define h s (t) = h(s,t). The homotopy is a homotopy of closed 
curves iff every h s , < s < 1, is a closed curve. A homotopy is continuously differentiable 
iff the map h is continuously differentiable. 

Two closed curves C , C\ in Q are homotopic iff there exists a continuously differentiable 
homotopy of closed curves h so that Cj is the image of hj, j = 0, 1. The domain Q is simply 
connected if every closed curve in Q is homotopic to the constant curve ("a point"). 

The intuitive meaning of "simply connected" is that every closed curve in Q can be shrunk 
to a point without leaving f2. Euclidean 3-space R 3 is simply connected, as is Euclidean 
3-space with the origin deleted. However Euclidean 3-space with the x-axis deleted is not 
simply connected. Another example: a cylindrical tank is simply connected, but a tank with 
a cylindrical hole up its axis is not, nor are donuts or coffee cups with handles. 
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Poincare's Lemma: Suppose that C R 3 is simply connected, and that that u is a 
continuously differentiable vector field on Q. Then the following hypotheses are equivalent: 

• u has zero circulation around any closed curve in Q; 

• V x u = 0; 

• u = V0 for a twice continuously differentiable scalar valued function <fi on Q. 

Assume that the volume occupied by a steady irrotational flow is simply connected. Then 
Poincare's Lemma implies that the velocity is the gradient of a scalar function, called the 
flow potential 0: 

V = V0 

This identity combined with the momentum balance law leads to Bernoulli 's Theorem, which 
you will develop in Problem Set 5. Expressing the incompressibility condition in terms of 
the flow potential leads to the Laplace Equation: 

3 d 2 S 

V 2 = V-V0 = £i4 = O 

i=l UX i 

arguably the most important equation in mathematical physics. 

To gain some idea of the powerful constraint the Laplace equation places on flow, suppose 
steady, irrotational flow takes place in a simply connected domain Q, and that the boundary 
confines the fluid. Necessarily, the fluid velocity is tangent to the boundary. That is, 

V0(x) • n(x) = v(x) • n(x) =0, xeffi 

This equation is properly known as the "no flow" boundary condition, or the (homogeneous) 
Neumann boundary condition. The pair consisting of the Laplace equation and the Neumann 
boundary condition is known as the Neumann Problem for 0. Use vector calculus identities 
from old homeworks and the divergence theorem to see that 

= / 0V 2 0= / (V • (0V0) - V0 • V0) 

= f 0V0-n- / |V0| 2 
Jan Jn 

= -[ |V0| 2 = -/ |v| 2 
Jn Jn 

The velocity field is continuously differentiable by standing assumption, so its length squared 
is continuous and nonnegative. Yet another important theorem from advanced calculus: if 
a nonnegative continuous function integrates to zero, then it is zero. 
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So we have shown: if steady irrotational flow is confined to a simply connected domain, 
then there is no flow at all!!! The velocity field vanishes identically, from which it follows 
that the external force field must be conservative (see Problem 2 below) and balance the 
pressure gradient precisely. The flow potential is constant - one cannot say more, since all 
of the conditions constraining it involve its derivatives. 

If you want steady nonzero confined flow, you will have to settle for nonzero circulation, 
or make it happen in a domain that's not simply connected! 

In fact we have proven more: if two steady irrotational flows in the same simply connected 
domain, with velocity fields Vi and v 2 , have identical normal components of velocity on the 
boundary, then the two flows have the same velocity field throughout the domain, and their 
pressure fields differ by a constant. Indeed 5v = v 2 — vi has zero curl and divergence also, so 
it is the gradient of a potential which solves the homogeneous Neumann problem. We have 
just seen that all such (twice continuously differentiable) potentials are constant, therefore 
the gradient 5v vanishes, i.e. Vi = v 2 . 

You can phrase this observation entirely in terms of potentials: if <f>i,<f>2 both solve the 
same inhomogeneous Neumann problem 

V 2 0i = / in tt, V0 • n = g on dtt, i = 1, 2 

defined by suitable right hand sides / : Q — > R and g : Q — >■ R, then 0i and 2 differ by a 
constant. That is, solutions of the Neumann problem are as unique as they can be. 

Another consequence: for steady irrotational flow in a pipe, as discussed earlier, the 
normal components of velocity at the inlet and outlet completely determine the flow at 
every point in the pipe. 

For very simple situations, this fact is enough to pin down the flow explicitly. Suppose 
for example that the pipe is straight, say with its axis parallel to the x\ axis, and simply 
connected. The inlet and outlet have the same area, and the normal (xi) component v of 
velocity is constant across the inlet and outlet and the same (as we saw, this last is required 
by conservation of mass). Then v(x) = (u, 0, 0) T is the gradient of 0(x) = uxi, which 
solves the Laplace equation and matches all of the boundary conditions. Therefore the flow 
potential is 0+ const., hence the constant v must actually be the velocity field in the pipe. 

Functions satisfying the Laplace equation are also called harmonic. Harmonic functions 
of two variables (say, x\ and x 2 ) are quite special: they are real parts of complex analytic 
functions of the complex variable x\ + ixi- Such potentials give rise to flow parallel to the 
(xi,x 2 ) plane (so called planar flow), and in special cases they can be determined explicitly 
by the method of conformal mapping, for which see practically any text on complex variable 
theory or elementary fluid mechanics. While few problems admit the conformal mapping 
approach, some that do are quite important. 
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4.5 Green's Functions 



This section shows how we may construct solutions u of the Poisson Problem 

V 2 u = <f) 

for a rather tame type of right hand side (f>. This is the beginning of a line of reasoning that 
leads eventually to construction of potential flows in various circumstances; here we merely 
sketch the beginning of these developments. 

First we need a few notions related to functions, and things more general than functions, 
with singularities. 

A test function on Euclidean <i-space H d is an infinitely often differentiable scalar- valued 
function which vanishes for large enough ||x||. The set of test functions is conventially 
denoted V; it forms a vector space. 

Note that any (mixed) partial derivative of any order, or linear combination of partial 
derivatives, or product with an infinitely differentiable function, of a test function is a test 
function. 

It's conventional to take for the scalar field of values the complex numbers, rather than 
the real numbers. All of our examples will be real-valued, but it makes very little difference 
to the theory developed below if complex-valued test functions are used. 

There is a natural notion of convergence of sequences of test functions: if {0 n }^L o is such 
a sequence, and denotes the zero function (which plays the role of the zero vector in the 
vector space V, then n — > if and only if (i) for some R > 0, 0„(x) = if |x|| > R, any n, 
and (ii) 

Qh+--+3d, 



SU P||x||<i? 

for each R > 0. 



dxl 1 ...dx J / 



-> 



It gets tedious to write the partial derivatives symbols in this expression, which motivates 
the definition: for a vector j = (ji, ...,jd) of non-negative indices, 

p)h+-+3d 



dx{ x ...dx J d d 



The set of non-negative integers is denoted Z+, so j G Z+. We can restate the second 
condition for convergence 4> n — > in V as: for any j G Zl, 



su P\\x\\<R 
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A distribution is a continuous scalar-valued linear function on the vector space V. The 
distributions also form a vector space in a natural way, denoted by V . 

The value of a distribution u E V on a test function G £> is conventionally denoted 
(u,4>). 

Continuity means that if u E V and {0 n }^L o is a sequence of test functions, then 

n (u, n ) ->■ 

Two examples of distributions are already familiar, if not under that name: 



The Dirac Delta function: Usually denoted 5, defined by the rule 

(5, 0) = 0(0) 

Caution: in many places the value of the delta function on a test function will be written as 
an integral: 

J rfx5(x)0(x) =0(0) 

Since the delta function is not a function, the integral notation for its value must actually 
be poetry - i.e. the left hand side of the preceding equation has no meaning, apart from 
denoting the right hand side. 



Integrable Functions An integrable scalar valued function g defines a distribution in a 
natural way, and it's conventional to use the same letter to denote the function and the 
distribution it defines: 

(9,4>) = J dxg(x)(f)(x) 

The overbar denoted complex conjugation, which of course is a no-op if g is real valued - so 
in that case I will drop the overbar. 

We extend various common operations on functions to distributions by transposition. 
Suppose that T : T> — > T> is linear and continuous. Just as for scalar-valued functions, 
continuity means that if n — > in the sense defined above, then T0 ra — > also. Define the 
transpose T' : T>' — > T>' by the rule: if u E V, E V, then 

(T' M ,0> = ( M ,T0> 



Differentiation, being continuous on T> (verify this!), extends to distributions via trans- 
position. That is if u E V, then u is the distribution defined by 



'AX A = ( w 

dx i) ' / \ dxi 
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This notion of derivative is often called the weak or distribution derivative. Note that if g 
and dg/dxi are integrable functions, then integration by parts shows that the (weak) (x, 
partial) derivative of g is actually the distribution defined by —dg/dxi. Thus the weak and 
usual notions of derivative coincide when both make sense. 

In view of this fact, it's conventional to confuse notations, and write — instead of the 

more cumbersome fgf^) u even when u is a distribution. Carrying on in this fashion, it's 
conventional to confused the transposed Laplacian on distributions with the usual Laplacian 
on functions, as the two coincide on distributions which happen to be (twice-differentiable) 
functions. 

A Green's Function for the Laplace operator is a distribution g e V for which 

V 2 g = 5 

(in which the left hand side is interpreted in the sense of weak derivative, of course). Clearly 
any two Green's functions, say g 1 and g 2 , differ by a solution of 

V 2 (3!-3 2 )=0 

i.e. g\ - ^2 is a "harmonic distribution". As it turns out, any harmonic distribution is 
actually (defined by) a harmonic function which is infinitely differentiable. We are about to 
see that the Green's functions themselves are definitely not so regular, as you would guess 
by looking at the right hand side of the defining PDE. 



Claim: g : R 3 ->■ R defined by 

<?(*) = ' 



47r||x|| 

is a Green's function for the Laplace operator. [Note: you should be able to call on your 
advanced calculus to show that the right hand side of this equation defines an integrable 
function (its singularity at zero is integrable), so this distribution is one of those defined by 
an integrable function.] 

To see this, we need to show that the defining equation is satisfied. Since both sides are 
distributions, this means showing that for any test function (f> e T>, 

(V 2 ^,0) = 0(O) 

This means (using the weak derivative definition several times) that 

<3,V 2 0) = 0(O) 

Gb£((x)V 2 0(x) 



/ 
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since g is actually an integrable function (hence its product with V 2 is integrable). 



The integral on the right hand side of the preceding equation is of course defined as an 
improper integral, since the integrand is not continuous - i.e. as the limit of proper integrals 

= rim e _> / cfe5f(x)V 2 (/>(x) 

■ / l|x||>e 

In this "punctured" domain of integration, the factors in the integrand are twice (actu- 
ally, infinitely often) continuously differentiable, so you can integrate by parts (i.e. use the 
divergence theorem: 

= lim e _> / dx V • (#V0) - Wg • V0 

J||x||>e 

= lim e _K, \ I dA(x) gV<f) ■ n - [ dx V • {Vgcfy + f dx V 2 g<p \ 

I J||x||=e ./||x||>e •'ll x ll> e J 

= lim e ^ \ [ dA(x) gV(f) • n - / dA(x) Vg ■ n0 + / dx V 2 g<p \ 

[7||x||=e J\\x-\\= e "'ll x ll> e J 

Now we can start getting rid of stuff. According to problem 3-5 (c), the last term on the 
right hand side drops out. Note that 



and that 



l x ll = e =>• 9(x) = 



[ dA(x) V0 • n 

J llxll=e 
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0(e 2 ) 



i.e. this integral is bounded by a multiple of the area of the sphere of radius epsilon, since 
the integrand is bounded. Thus 



dA(x) gV(p ■ n 



= 0{e) 



so the first term goes to as e — > 0. 
Finally, note that 

• n(x) = — x/||x|| (this is the outward unit normal of the outside of the ball of radius e!); 

• from problem 3-5 (a) and a little algebra, 

Vs(x) = 



so 



V#(x) • n(x) 
on the boundary ||x|| = e.; 



47r||x|| 3 
1 



1 



47r||x|| 2 Aire 2 
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• if ||x|| = e, then 0(x) = 0(0) + 0(e).. 

Putting all of this together, 

-/ dA(x)Vg -n(f) = -(0(0) + 0(e)) / dA(x) Vg ■ n 

J||x||=e ./||x||=e 

= 0(0) + 0(e) 0(0) 

as e — )• 0. q.e.d. 

The Green's function bears the same relation to the Laplace operator that the inverse of a 
matrix bears to the matrix. To see this, we need some more information about transposition. 

The translation operator T y : V — > V moves its inputs by y G R n : 

(T y 0)(x) = 0(x + y) 



Claim: for any G V, y G R n , i — 1, n, and sequence of reals h k — > 0, 

in the sense of convergence in T>. 

Proof: we must show that (a) all of the functions on the left of the above limit statement 
vanish outside some large ball, and (b) they and their derivatives of all orders converge 
uniformly to the RHS. Concerning (a), suppose that vanishes outside a ball of radius R, 
and \hk\ < H for all k (there must be some such H, since the sequence is convergent). Then 
the function on the left hand side of the above equation vanishes outside the ball of radius 
R + H+\y\. 

According to the Mean Value Theorem, for each k there is h k between h k and so that 
for any x G R n , 



l(T y+ ^ ei 0-T y 0)(x)-T y ^(x; 



d(j> 



0(x + y + h k ei) - 0(x + y) 



Since the is uniformly continuous (being continuous and vanishing outside a bounded set), 
the last absolute value tends to zero uniformly in x as k — > oo. 

If G V and u G V, then the rule 

y^( M ,T y 0> = (7>,0> 



50 



defines a function on R d , and that this function is infinitely often differentiate (proof: use 
continuity and the claim to see that any first partial exists; now apply the same reasoning 
repeatedly, using the fact that any mixed partial derivative of a test function is a test 
function.) This function has a name: it is the cross correlation of u and 0, written 

u*(j)(y) = (u,T y( f)). 

As a by-product of the induction, you get the rule: for any j G Z|, 

D*u * <fi — u * D*<j) = (jy)'u -k <f) 



Now take for u the Green's function g derived above, a distribution defined by an inte- 
grable function. Applying the second partials in turn and adding them up (distributions are 
linear functions!), we get 

V 2 s * 0(y) = (V 2 )'<7 * 0(y) = 0(y) 
= (5,T y 0) = (T y 0)(O) = 0(y) 

for any y G R d , i.e. 

V 2 ((?*0) = 0. 

That is, the crosscorrelation of the Green's function with the right hand side (assumed to 
be a test function) solves the Poisson equation! 

Recall that the Green's function is actually a function, i.e. acts as a distribution via 
integration. Also, it is an even function of its argument. So 

9 * 0(y) = J dxg(x)<f>(x + y) = J dxg(x - y)0(x) = J dxg(y - x)0(x) 

The last integral defines another common type of product of functions, namely the convo- 
lution: the convolution of an integrable function / and a test function ip e V is the function 
f * 4> defined by 

y ^ / * 0(y) = j dxf(y - x)0(x) 
It's easy to verify that for integrable /, ip e V, 

f*<P = /*0= /*</> 

where <f>\-> <f> is the parity operator (thanks to Rani Namour for suggesting this name): 

0(x)=0(-x). 

The definition of the parity operator makes sense for any type of function - test, integrable, 
or neither. 
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This observation also motivates the extension of the notion of convolution to distributions: 
by definition, if u G V, G V, then u * is the (C°°) function given by 

u * 0(y) =< TyU, <j)>=u-k 0(y). 

For the solution of the Poisson problem with right-hand side G V, we can just as well 
write 

V 2 {g * 0) = V 2 {g * 0) = V 2 {g * 0) = V 2 {g * 0) 

since g is even (i.e. g = g). In other words, you get the solution of the Poisson problem by 
convolving the right hand side with the Green's function. 

You should be aware of a common notational convention in the applied literature. The 
application of the transpose translation to the delta distribution produces another distribu- 
tion, and just as with delta it's common to denote it as if it were a function: 

5(x-y)" = w (7j*)(x) 

which simply means that 

(<H--y),0> = 0(y) " = / ^5(x-y)0(x)" 

As before, the last "equation", giving a common notation for the evaluation of this distribu- 
tion, has no independent meaning. 

It's common to use this notation to express the solution of the Poisson problem by the 
Green's function in the form 

V 2 (?(x-y) = 5(x-y). 



4.6 Problems 



Problem 1: Suppose that M is a continuously differentiable matrix- valued function on a 
subset of R 3 . Show that the divergence V • M is the only vector valued function defined on 
the same domain as M and satisfying 

(V • M) ■ c = V • (M T c) 

for an arbitrary constant vector c G R 3 . 
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Problem 4-2: Suppose that v(x, t) is the velocity field of an ideal fluid moving under the 
influence of a conservative external force density f , i.e. 

f(x,t)=p V/3(x,t) 

for some scalar field j3 (which, you may note, must have units of force x distance per unit 
mass, or energy per unit mass). Suppose that a material part of the fluid occupies the volume 
Q t at time t, and that the velocity field is tangent to the boundary at all times and places: 

v(x,t) • n(x,f) = 

for any point x on the boundary of Qt with normal n(x, t). 

Show that the kinentic energy of this part is constant in t, i.e. that 

at 2 Jn t 

Hint: You will have to remember the transport theorem, and work out for example that 
V|v| 2 = 2(Vv) T v. 



Problem 4-3: Prove the following special case of Bernoulli's theorem: suppose that 
v(x),p(x) are the velocity and pressure filed of a steady and irrotational motion of an ideal 
fluid moving under the influence of a conservative external force f = poVj3. Show that 




is constant, ie. independent of position. 



Problem 4-4: Suppose that v(x),p(x) are the velocity and pressure field of a steady and 
irrotational motion of an ideal fluid moving under the influence of no external force, i.e. in 
my usual notation f = 0. Suppose that the volume occupied by the fluid excludes a domain 
E which is filled with a hard, nonfluid obstacle. Show that the total force exerted on the 
obstacle by the fluid is exactly 



Po r 

2 J a 



tL4(x)|v(x)| 2 n(x) 

an 



Hint: You will need to use the result of Problem 3, and also justify several bits of calculus, 
for instance the identity 

f dA(x)n(x) = 

for any bounded domain Q - this is a consequence of the divergence theorem for matrix- valued 
functions. 
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Problem 4-5: Show that 

<?(x) = ^log(||x||), x G R 2 
is a Green's functions for the Laplacian in R 2 . 
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Chapter 5 

Elastic Fluids and Acoustic Waves 



5.1 Characterization of Elastic Fluids 



While everyday fluids are mostly very difficult to compress, hence well modeled by incom- 
pressible continua, gases are notably compressible. In another respect, gases behave like 
(ideal) fluids: they support only small shear stress. In the idealized limit of no shear stress, 
gases might be modeled as compressible, or elastic, fluids: once again, the stress is scalar 
(pressure) but now motions in general do not preserve volume, so the divergence of the 
velocity field is no longer constrained to vanish. 

This model consists in part of the laws of mass and momentum balance: 

§ + (V.v)p = 

These conditions are not sufficient to determine the motion uniquely under any reasonable 
circumstances - a replacement is needed for the velocity divergence constraint. The replace- 
ment is the constitutive law relating density and pressure: 

V = p{p) 

This law is supposed to reflect on a macroscopic scale the molecular properties of the material. 
For real gases, it is observationally true that an increase in density is correlated with an 
increase in pressure, so assume in addition that 

dp 
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In an ideal fluid, it was not possible to compute the pressure from the density, as the 
density was assumed constant but the pressure was free to vary from place to place (in 
general). Moreover the pressure did not determine the velocity field either, as one could add 
a constant of arbitrary size to the pressure without violating any of the equations of motion. 
This freedom no longer exists for an elastic fluid. 

A familiar example of an elastic fluid constitutive law is the ideal gas law, which relates 
pressure p, density p, temperature 9, and a constant R characteristic of the material: 

p = RQp 

If you assume that a volume V is filled with a material obeying the ideal gas law, with 
homogeneous density, temperature, and pressure, then 

pV = R6pV = RmO = const 9 

which is the famous "PV=T" of high school chemistry, apart from choice of units. 

There is a Bernoulli Theorem for elastic fluids, just as for ideal fluids - see [Gurtin, 1981], 
pp. 130 ff. I want to emphasize two other aspects of typical gas behaviour here, both of 
which will be developed much more extensively in the second half of this course. 



5.2 Acoustics 



If you set p(x, t) = po, a constant, v(x, t) = 0, and p(x, t) = po = P{po), and assume absence 
of external force, you have a constant or equilibrium solution of the elastic fluid equations of 
motion: pressure and density are homogeneous and the material (gas) is not moving. Now 
suppose that you add a very small, possibly time-varying force density <5f to the right hand 
side of the momentum balance law. Then presumably the motion varies by a small amount 
from the equilibrium solution: 

V = Po + Sp, p = p + 5p, v = 5v 

where anything with a 5 in front of it is supposed to be small, and even to have small 
derivatives. 

Substitute these expressions in the mass and momentum balance laws and the constitutive 
law: 

*e£M + (V.*r)<A, + «,) = <> 
dSv 

Po + Sp)-zr + (V5v)5v + V(p + Sp) = Si 
at 

p + 5p = P(p + Sp) 



56 



Presumably when all things with 5s are small, the products of things with 5s are much 
smaller. So approximate the balance laws by dropping these product terms, and the consti- 
tutive law by using its first order Taylor series. Make use of the equilibrium equations to 
eliminate the terms without 5s. The result is the system 

^ + p V • 5v = 

p0 ~df + p = 

Since dP/dp > 0, you can solve this equation for 5p in terms of 5p, and substitute the result 
in the conservation of mass. The rate of change of pressure with respect to density is so 
important that it gets its own name: 



IdP . , 



is the sound speed in the fluid (the reason for the name will become apparent shortly). 

With this substitution to eliminate 5p, the remaining two equations become the acoustic 
equations of motion: 

1 05p _. 
p Q c z at 

Po ~cW + P = 

This system is equivalent to a single second order partial differential equation: differen- 
tiate the first equation with respect to time, use the second equation to eliminate the time 
derivative of velocity, and you get 

1 d 2 5p „ ( 1 r \ V-5f 
v5p 



p c 2 dt 2 \po J p 

Since p is constant, this is equivalent to 

1 d 2 5p 



c 2 dt 2 

which is the acoustic wave equation. 



V 2 5p = -V • 5f 



5.3 Problems. 
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Problem 5-1: [50 pts] This problem concerns small amplitude motion of an elastic fluid 
contained in a domain Q C R 3 . The motion is described by the perturbational fields 
5p(x, t), 5v(x, t). These are assumed to satisfy the acoustic wave equation described in 
section 5.2. Define the function of time 

|2 <5p(x,t) 2 l 



Show that £"(£) is independent oft if no external force (perturbation) is acting on the elastic 
fluid and either 

(i) <Sv(x, t) -n(x) = 0, x G dQ (the "no-flow" condition - i.e. Q is enclosed by an impermeable 
barrier), or 

(ii) 5p(x.,t) = 0, x G d£l (the "free surface" or "pressure release" condition). 

As will be established in Ch. 7, E(t) is actually the energy in the perturbational field, and 
this result can be paraphrased by the statement that this energy like quantity is conserved 
by the motion, provided that either no-flow or pressure-release conditions are imposed on 

on. 



Problem 5-2: [25 pts] Show that any function of the form 

6p(x,t) = f(t - p -x) 

is a solution of the acoustic wave equation, for any waveform f and any slowness vector p, 
provided that 

, , 1 
P = — 

c 

(the units implied by this equation explain the name of p). 



Problem 5-3: [25 pts] Suppose that f2 = {x : rr 3 > 0} is the positive-a^ half space, and 
that f+(t) and p + are a waveform and slowness vector respectively. Exhibit second pair 
f-(t), P- of waveform and slowness vector so that the sum of the corresponding solutions of 
the acoustic wave equation 

<5p(x, t) =/+(*- p+ • x) + f-(t ~ P- • x) 

is also a solution of the acoustic wave equation, for which the free surface condition holds: 

Sp(x, t) = 0, x G dQ 

["Exhibit" here means: give explicit formulas for /_,p_ in terms of / + ,p + .] 

Show that the vectors p + , p_ make equal angles with the boundary normal e 3 . [This result 
is known as SnelVs law of reflection] 
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Chapter 6 



Newtonian Fluids and the 
Navier-Stokes Equations 



6.1 Basic Consitutive Assumptions 



A Newtonian fluid is a material whose motion satisfy the following constraints: 

• density is constant, p(x, t) = p ; 

• only volume-preserving motions occur, i.e. V ■ v = 0; 

• the stress is a pressure plus a linear function of the velocity gradient Vv, which measures 
the relative velocities of nearby fluid particles. 

Thus Newtonian fluids are incompressible and have homogeneous density, like ideal fluids. 
Unlike ideal fluids, Newtonian fluids can support shear stress, when the velocity is nonuni- 
form, i.e. changes in space. 

The stress is (in part) a matrix-valued linear function of a matrix, and such a thing is a 
bit complex. You can express it with the help of a four- index array, C^: 



3 



dv k 
dxi 



k,l=l 



or more compactly 



S = - 



P I + C[Vv] 
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Together with the laws of conservation of mass and momentum and the volume-preserving 
condition in their familiar forms, this condition provides the ten equations one would expect 
to need to determine the ten unknowns p, v, and S. 

Stress is symmetric, hence has six independent components; the velocity gradient appears 
to have nine, so there would appear to be no less than fifty-four independent constants Cijki 
involved in the relation between stress and velocity. That there are actually not nearly so 
many is a consequence of a natural physical invariance property, to which we now turn. 



6.2 Observer independence 

The principle of observer independence codifies the idea that the constitutive assumptions 
characterizing a type of material should not depend on how it is observed, i.e. on changes 
of coordinates, even time varying ones, which are distance-preserving. This is not to say 
that the various component functions which go into the statements of the dynamical laws 
(conservation principles) will necessarily remain the same, only that so long as distances 
remain the same in a new coordinate system, the form of these equations will remain the 
same as in the old coordinate system. 

In equations, imagine that you have two coordinate systems, x and x', related by a 
time- varying Euclidean transformation (which preserves distances!): 



in which Xq : R — > R 3 and R : R — > R 3x3 are differentiable and R(t) is orthogonal for 
each t and det R = 1. As pointed out in section 1.4.2, only afline coordinate changes with 
orthogonal R preserve all distances. You'll see why det R must be positive shortly. 

You are to think of x and x' as two different descriptions of the same point in space, one 
of whom is moving and spinning, possibly with time-dependent rates, relative to the other. 

If x(X, t) is a motion, then the same motion, described in the x' coordinate system, is 



x' = x'(x,t) = x' (t) + R(t)(x - x ) 



x'(X, t) = x' (t) + i*(t)(x(X, t) - x ) = x'(x(X, t),t) 



Various quantities related to the motion functions are also related, as follows. 



From the chain rule, 



<9x' 
~dt 



(X,t) 




dt 



(t) + ^(t)x(X,t) + R(t)%(X,t) 



V x x'(X,t) = V x x'(x(X,t),t)V x x(X,t) = i*(t)V x x(X,t) 
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This last formula justifies the limitation det R > 0: otherwise the essential property det Vx > 
cannot hold in both coordinate systems. Write X'(x', t) for the inverse function of x'(X, t), 
i.e. X'(x'(X, t), t) = X for any X G R 3 . If x' and x refer to the same material, i.e. if 
x' = x'(x, t), then necessarily X'(x'(x, t), t) = X(x, t), since the material coordinates simply 
label the material. 

The material velocity field v'(x',t) in the x' coordinate system, is 

v'(x^ = ^(X'(x',0,0 

so 

<9x' 

v'(x'(x,t),t) = ^ r (X'(x'(x,t),t),t) 
= -^(Xfx, t),t) = - + -^(«)x(X(x, (),() + B(0 w (x(x, 0,0 

= #W + §W* + *«v<*,<) 

So the velocity in the x' coordinate system differs from that in the x coordinate system both 
by the instantaneous rotation R(t) and by the velocity of the x' coordinates relative to the 
x coordinates (the first two summands). 

The velocity gradients (Jacobians) are also related: 

V x v'(x'(x,t),t) = (V x ,v')(x'(x,t),t)V x x(x,t) 

= (V x V)(x / (x,t),t) J R(t) 
dR 

= i2(t)V x v(x, *) + —(*) 

or 

(V x ,v')(x'(x,t),t) = R(t)VM^t)R(t) T + ^(t)R(t) T 
which is the required relation. 

To see how velocity divergences are related, note the following facts: 

• a matrix and its transpose have the same trace; 

• skew symmetric matrices Q (Q T = —Q) have zero trace: trQ = trQ T = — trQ, so 
trQ = 0; 

• the trace of a matrix product does not depend on order: tiAB = trBA so long as both 
sides make sense and are square; 

• similar matrices have the same trace: tr(RAR~ r ) = ti^AR^R) = trA 
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if R(t) is differentiable and orthogonal for each t, then ^(t)R T (t) is skew-symmetric: 



The velocity divergence is the trace of its Jacobian: 

(V x / •v / )(x / (x,t),t) =tr(V x V)(x'(x,t),t) = 

= tr J R(t)V x v(x,t) J R(t) T + tr (^(t)R(t) T ) = V x -v(x,t) 



using the just-listed results on traces. 



One prominent quantity has not been mentioned so far in this discussion: the stress. We 
have actually already seen how the stress changes when coordinates change by a Euclidean 
transformation (section 1.4.2), and that analysis applies equally well when the Euclidean 
transformation depends on time: 

S'(x!(x,t),t) = R(t)S(x,t)R(t) T 



These relations quickly imply that the constitutive assumptions appearing in the fluid 
models of the last two sections are independent of observer (homework!). For Newtonian 
fluids, two critical properties of the mapping C follow from observer independence: 

1. C depends only on the symmetric part of the velocity gradient: for any 3x3 matrix 

C[L] = C 



2. C is an isotropic transformation, i.e. for any symmetric L, orthogonal R, 

RC[L]R T = C[RLR T ] 



Proof of (1) : First, observe that for any 3x3 matrix L, there is a motion x(x ,t) for 
which v(x, t) = Lx, hence Vv(x, t) = L: set 



x(X,t) = e* L X. 
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Note that a(t) = det e tL satisfies 

= (trL)a(f) 

and a(0) = 1, whence det e tL is always positive. The map X h-> x(X, t) is plainly invertible, 
so in fact x(X, i) actually defines a motion. 

Note that 

^(X,t) = Le' L X = Lx(X,t) 

whence v(x, t) = Lx as claimed. 

Note that if x(X, i) as just described is to be the motion of a Newtonian fluid, it must be 
volume-preserving, whence V ■ v = trL = 0. So assume from now on that L has zero trace. 

The stress tensor for a Newtonian fluid under this motion must take the form 

S(x,t)=p(x,t)I + C[L] 

where p(x, t) is the pressure, which we leave unspecified, and tr C[L] = 0. 

For any 3x3 matrix L, motion x(X, t) as just constructed, denote by W the skew- 
symmetric part of L: 

W= 1 -(L-LT) 

Set 

R(t) = e~ tw , x'(x,t) =R(t)x 
Since W is skew-symmetric, R(t) is orthogonal for each t. 

In the x' coordinate system, 

1 S / (x / ,t)=p / (x',t)/ + C[V x V(x',t)] 

= R(t)S(x, t)R(t) T = R(t)(p(x, t)I + C[L])R{t) T 

= p(x,t)I + R(t)C[L]R(t) T 

Since S and S' are similar, they have the same trace, so p'(x',t) = p(x, t). As shown 
previously, 

V x ,v'(x',t) = J R(t)V x v(x,t) J R(t) T + ^r(t)R(tf 
= R{t)LR{t) T - W 
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At t = 0, then, 



S"(x', 0) - p'(x', 0)/ = C[L - = C 
= S(x,0) -p(x,0)7 = C[L] 



q.e.d. 



Proof of (2) : this time take R constant, x' = Rx, so that 

V x V(x,£) = RV x v(x,t)R T 

From (1), we might as well confine our attention to symmetric, zero-trace matrices L: any 
matrix orthogonally symmetric to such a matrix is also symmetric and zero-trace. For the 
motion generated by such a matrix as above, and arbitrary orthogonal R, 

S"(x', f) - p'(x', t)I = C[RLR T ] 

= R(S(x,t) -p{x,t)I)R T = RC[L]R T 

q.e.d. 

What does isotropy mean? Recall that C[Vv] is (part of the) stress: it transforms under 
a coordinate change x' = x' + Rx. (with orthogonal R) like any other stress, namely to 
RC[Vv]R T . That is, if you measure this part of the stress tensor in the new coordinate 
system x' ; it will be related to the same part measured in the old system as indicated. Thus 
you can take your new measurements C and convert them back to the old measurements by 
the inverse similarity transform: C = R T C'R. On the other hand, the velocity Jacobian Vv 
also transforms by similarity, Vv' = RVvR. So the following two sequences of operations: 

• measure velocity Jacobian in coordinate system x 

• compute corresponding stress 

• transform stress to new coordinate system 

and 

• measure velocity Jacobian in coordinate system x 

• transform velocity gradient to new coordinate system x' 

• compute stress in new coordinate system 



will give you the same result. In that sense, orientation of the coordinate system does not 
make any difference in the prediction of stress from velocity Jacobian. 
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6.3 Isotropic Transformations 



The structure of isotropic transformations on the space of symmetric matrices is quite re- 
stricted. There is no reason to limit ourselves to three dimensions in studying the implications 
of isotropy, though the only applications in this course will be three dimensional! 

Structure Theorem for Isotropic Transformations: Suppose that C is a linear trans- 
formation on the vector space of symmetric n x n matrices, and is isotropic: for any orthog- 
onal n x n matrix P, 

C[RLR T ] = RC[L]R T 
for any symmetric L. Then there exist constants /i, A so that 

C[L] = 2/2L + A(trL)/ 
Proof: The spectral theorem for symmetric matrices states that 

k 

i=i 

where the Aj are the k (real) distinct eigenvalues of L and Pj corresponding orthogonal 
projectors onto the subspaces of eigenvectors for Aj. The dimension of range (Pj) is the 
dimension of the eigenspace for Aj, i.e. its multiplicity. If the multiplicity of Aj is 1, then 
Pj = uu T where u is a unit eigenvector with eigenvalue Aj. In general Pj is a sum of rank 1 
matrices of this form, with u ranging over an orthonormal set of eigenvectors for Aj. 

First, we claim that every eigenvector of L is an eigenvector of C[L\. Let u be an 
eigenvector of L, and choose R to be the reflection across the plane orthogonal to u, i.e. 

Ru = — u, Pv = v if v • u = 

Since R has the same eigenvectors as L, R and L commute, hence 

RC[L]R T = C[RLR T ] = C[LRR T ] = C[L] 

i.e. 

RC[L] =C[L]R 

Thus 

RC[L]u = C[L]Ru = -C[L]u 

which means that C[L]u is an eigenvector of R with eigenvalue —1 - which is only possible 
if C[L]u is a multiple of u, since the —1 eigenspace of R has dimension 1. In other words, u 
is an eigenvector of C[L}. Since u was any eigenvector of L, this establishes the claim. 
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Since C is linear, its behaviour is therefore determined by its action on matrices of 
the form L = uu T , where u is a unit vector. The eigenvalues of L are 1 (multiplicity 
1) (eigenvectors = multiples of u) and (multiplicity n — 1) (eigenvectors = all vectors 
orthogonal to u). Thus u is also an eigenvector of C[L], and the entire subspace of R n 
orthogonal to u also consists of eigenvectors. If all vectors in a subspace are eigenvectors, 
then they all belong to the same eigenvalue. So 

C[L] = a(u)uu T + (3(u)(I - uu T ) = 2//(u)uu T + A(u)J 

where the eigenvalues of C[L] are a(u) and /3(u), which so far as we know at the moment 
depend on u, and fx = |(a(u) — /3(u)), A(u) = /3(u). 

To see that in fact A and /i are independent of u, choose any orthogonal matrix R and 
set v = Ru. Then 

C[vv T ] = 2/x(v)vv T + A(v)J 
= C[Ru{Ruf] = C[Ruu T R T ] = RC[uu T ]R T 
= R(2fi(u)uu T + X(u)I)R T = 2/i(u)vv T + A(u)7 

so 

(2(/i(v)-Mu)))vv T = -(A(v)-A(u))/ 
Since vv T and I are linearly independent, it follows that 

/i(v) = /i(u), A(v) = A(u) 

As R ranges over all orthogonal matrices, v = Ru ranges over all unit vectors. Therefore in 
fact A and /i are independent of u. 

Since tr uu T = 1 for any unit vector u, we have actually shown that 

C[uu T ] = 2/mu T + A(truu T )J 

Since C is linear, any symmetric L is a linear combination of orthogonal projectors 
(spectral theorem!), and any orthogonal projector is a linear combination of rank 1 matrices 
of the form uu T , the "if" part of the theorem follows. 

The "only if" part is a homework excercise. q.e.d. 

Corollary: Special case of symmetric zero trace matrices: Suppose that C is an 
isotropic linear transformation on the vector space of symmetric zero trace matrices. Then 
for some constant /i, 

C[L] = 2/iL 



66 



Proof: Again, the proof of the "if" part is given here. The "only if" part is actually just 
a special case of the "only if" part of the previous theorem. 

For an arbitrary symmetric matrix L, define 



C[L] = C 



n 



(trL)J 



This makes sense, since the matrix in brackets on the right has trace zero. It's very easy to 
verify that C is isotropic, so by the preceding theorem 

C[L] = 2pL + \(tr L)I 

for suitable constants p, A. If L has zero trace, then C[L] = C[L], so the corollary is proved, 
q.e.d. 



6.4 The Navier-Stokes Equations 

Over the last two sections we've seen that the constitutive assumption regarding the form 
of stress in Newtonian fluids, 

S = -pI + C[Vv] 

together with the requirement of independence of observer, mandate that for some constant 
p, in fact 

S = -pI + /i(Vv + (Vvf) 

Note that V • (Vv) T = V(V • v) (exercise!) and the latter vanishes because all Newtonian 
fluid motions are volume preserving. Also V • (Vv) = V 2 v, where the right hand side 
means the vector of Laplacians of the velocity components. Since the const ant- density and 
volume-preserving assumptions imply the conservation of mass, the characteristic equations 
of Newtonian fluid motion are thus 

Po^ + Vp-/iV 2 v = f 

V-v = 

These are the Navier-Stokes equations. Given po, p, and f, they consitute four nonlinear 
partial differential equations for the four components of p, v. 

The scalar p is the coefficient of viscosity: it regulates the strength of shear (viscous) 
internal forces in the fluid. 

The number of equations appears to be adequate, so that the Navier-Stokes equations 
should determine a Newtonian fluid motion given its initial state. On the other hand, it is 
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far from obvious that these equations have solutions. It is worth emphasizing once again 
that one cannot appeal to physics to assert the existence of solutions: many (reasonable but 
somewhat arbitrary) approximations and assumptions have gone into the formulation of the 
Newtonian fluid model, so that it is clearly only an approximation of the actual physics of 
real fluids. So the fact that real fluids flow does not imply that the Navier-Stokes equations 
have solutions! These equations, like the other models presented in this chapter, are models, 
and nothing more: their existence is justified if they predict actual motions of real materials 
well under the limited range of circumstances for which they are intended. Existence of 
solutions is a precondition for making reasonable predictions, and must be established as 
a property of the model. Solutions of the Navier-Stokes equations exist for some classes 
of initial data and limited time intervals. Global existence of solutions is still a subject of 
research. 

If the motion is steady (all time derivatives are negligible) and slow (velocities are small, 
so quadratic quantities in the velocity, notably (Vv)v, are negligible), then the Navier-Stokes 
equations reduce to the Stokes equations: 

Vp - /iV 2 v = f 

V-v = 

Take the divergence of the first equation: since the divergence commutes with the Laplacian, 
all that's left in view of the second equation is 

V 2 p = V • f 

With f given, this is a Poisson problem for p, i.e. an inhomogeneous version of the Laplace 
equation. 

The second chapter of these notes shows how to solve the Poisson problem for various 
boundary conditions modeling natural constraints on the fluid. Assuming that we can solve 
it, we have p. Then the Stokes equations themselves might be viewed as three more Poisson 
problems for the three components Vi of v. Thus solving the Stokes equations amounts to 
solving four Poisson problems. 

[Gurtin, 1981], sections 22-24, describes some basic properties of solutions of the Navier- 
Stokes equations. 



6.5 Problems. 



Problem 6-1: Show that the ideal fluid, elastic fluid, and Newtonian fluid models satisfy 
observer independence. What this means: you need to show that the constitutive assumptions 
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of each model are preserved under "change of observer", i.e. time-dependent Euclidean 
transformations of the form presented in the text: x'(x, t) = Xq(£) + i?(i)(x — x ) with R(t) 
orthogonal for each t, det R(t) = 1, and Xg,i? continuously differentiable as functions of t. 
The constitutive assumptions are: 

• ideal and Newtonian fluids: the density is constant 

• ideal and Newtonian fluids: the motion is volume-preserving 

• ideal and elastic fluids: the stress tensor is scalar (S = —pi) 

• elastic fluids: the pressure is a monotone increasing function of the density 

• Newtonian fluids: S = —pi + /i(Vv + (Vv) T ) 

You need to show that each of these assumptions is also true in the x' coordinate system if it 
is true in the x coordinate system. For the assumptions about the form of stress, remember 
that the desired result is 

S'(x!(x,t),t) = R(t)S(x,t)R(t) T 

In showing independence of observer for the Newtonian fluid model, you will incidently show 
the "only if" part of the theorem on the structure of isotropic transformations. 



Problem 6-2: Once again suppose that coordinate systems x and x' are related by a 
Euclidean transformation: x'(x, t) = x(,(t) + i?(t)(x — x ), with x(,(t) and R(t) as before. 
Show that the conservation of momentum law in the x' coordinate system takes the same 
form 

DV 

as in the x coordinate system: 




= V x • S + f 



and state an explicit relation between the external force vector fields f and f . Observe that 
f (x'(x, t),t) = i?(t)f(x, t) (as one might expect from the covariance of force) only if 
are independent of t. Changes of observer satisfying this constraint are called Galilean. 



Problem 6-3: Assume that po,v are a constant scalar resp. vector. Observe that they 
represent a slow, steady motion of a Newtonian fluid, i.e. solve the Stokes equation with van- 
ishing external force (f = 0). Suppose that p(x, t) = p + Sp(x, t), v(x, t) = v + 5v(x, t), f = 
5f represents a small perturbation of this constant solution. Following a procedure similar to 
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that used to obtain the acoustic equations from the elastic fluid model, derive the linearized 
Navier-Stokes equations under the assumption that all 5 quantities are small: 

Po~gf + Po(V5v)v + V5p - /iV 2 5v = 6f 
V • 5v = 
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Chapter 7 
Linear Elasticity 



Stress in elastic materials is a function of the Jacobian of the motion function: 

S(x,i)=«S((V q x)(q(x,f),f)) 
The function S : R 3x3 — > is characteristic of the elastic material. 

The meaning of this assumption is partly in what it rules out: in contrast to Newtonian 
fluids, elastic material stress does not depend on the relative velocities of nearby particles, 
but only on their relative displacement. 



7.1 Independence of Observer 

The change of variables 

x'(x,t) =x' (t) + i?(t)x 
with R(t) orthogonal, as in the last Chapter, gives through the chain rule 

V q (x'(x(q, *),*)) =i?(t)V q x(q,t) 

So the independence of observer property of stress, also explained in the last Chapter, be- 
comes 

S'(x!(x,t),t) = R(t)S(x,t)R(t) T = i?(t)5((V q x)(q(x,t),t)) J R(t) T 
= S(i2(f)(V q x)(q(x,t),t) 
Since Vx can take any value in R 3x3 , this condition implies that 

RS(L)R T = S(RL) 
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This condition is equivalent to the independence of observer property. 

One immediate consequence is that the stress is determined by the symmetric part of the 
Jacobian: since any matrix L has a polar decomposition L = RQ with orthogonal R and 
symmetric Q, 

S(L) = S{RQ) = RS(Q)R T 
Thus if you know S on Rg*m m , you know S on all matrices. 

Irritating Question: The independence of observer property mandated that deti? = 1; 
however that property is not assured of the orthogonal factor of the polar decomposition. 
Why does this fact not invalidate the claim just made (that S is determined by its action on 
symmetric matrices)? (Hint: what is the relation between det R and det(— R) if R is 3 x 3 
(or odd x odd)? 



7.2 Isotropy 

An elastic material is isotropic if the stress produced by a motion remains the same if the 
material is rotated before the motion. That is, 

5(V q x(q,i))=5(V q x(q',*)) 

where q' = qo + -R T (q — qo) are the coordinates of the point at q in the original reference 
configuration after it is rotated about qo by the rotation matrix R T . Note the ordering on 
the right hand side: points q in the reference configuration are first rotated (q — > q'), then 
moved (q' — > x) with the same motion function x. [Using R T rather than R to denote the 
arbitrary rotation matrix figuring in this concept is just a notational convenience, as you 
will see.] According to the chain rule, V q x(q',t) = V q 'x(q', £)V q q' = V q /x(q', t)R T . Set 
L = Vx(q ,t); since q' = qo when q = q , the equation above reads 

S(L) = S(LR T ) 

Since L may be any matrix (why?), this equation characterizes the stress function of isotropic 
elastic materials. 

Note that independence of observer also implies that 

RS(L)R T = S(RL) = S(RLR T ) 

Slightly generalizing the definition of isotropic linear transformation given in the last Chap- 
ter, then, you could say that isotropic materials have isotropic stress functions. 
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Some materials behave in a nearly isotropic elastic fashion, at least for small deformations 
- some plastics, some metals, sandstones, unreinforced concrete, amorphous salt. Others do 
not, exhibiting considerable change in stress after rotation of the reference configuration - 
muscle, composites (fiberglass, carbon fiber reinforced epoxy, reinforced concrete), shales. 



As was the case for elastic fluids, the theory of elasticity simplifies considerably when motions 
are small enough to neglect higher order effects. For variety, this section carries out the 
linearization process in a slightly more formal way than did Chapter 5. 

The easiest way to express the small- motion hypothesis is in terms of a 1 -parameter family 
of motions X/ l (q, t) = q + /m(q, £). Here h is a small parameter, and u is the infinitesimal 
deformation. We will seek the dominant parts of all quantities by expanding them in powers 
of h and dropping terms of order higher than linear in h. These expansions are legitimate 
when the quantities being expanded have enough derivatives, i.e. are well approximated by 
their Taylor series. These notes make the standing assumption (throughout this chapter 
only!) that plenty of derivatives are always available. 

To begin with, 



where the ellipses denote the remainder in this first order Taylor series, ie. "quantities of 
second order and higher" in h. 

S(I) is the stress at zero displacement, called the residual stress in the material. 

For simplicity, these notes assume that the residual stress vanishes. Not all elastic mate- 
rials have vanishing residual stress: as the name suggests, prestressed concrete is an example 
of material deliberately configured to have nonvanishing S(I). 

In coordinates, the linear (order 1 in h) part of the stress is 



where the four-index array C^i is the Hooke tensor. 

Since the stress must be symmetric (so that angular momentum balance is maintained), 
clearly C ijkl = C jm , all k, I. 



7.3 Linear Elasticity 



S{Vx h ) = S(I + Wu) 
-S(I) + h—(I)Vu + ... 
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For brevity, we will write 

so C is a linear transformation: R 3x3 — > R, 3 y^ m , with "matrix" Cijki- 

Actually the linear part of the stress depends only on the symmetric part e of Vu. To 
see this, note that for any matrix L, 

e hL = I + hL + ... 

Write Vu = e + r, where e = |(Vu + (Vw) T ),r = ^(Vw — (Vu) T ). Then 

S(Vx h ) = S(J + /i(e + r)) 

= S((I + hr) + he) = S((I + hr)(I + he) - h 2 re) 
= S(R(I + he) + ...) 

where R = e hr is a rotation matrix because r is skew-symmetric. It is easy to see that the 
higher order terms (...) can only contribute the higher order terms in the Taylor expansion 
of S: 

= S(R(I + he)) + ... = RS(I + ht)R T 
(because of observer-independence) 

= (/ + hr + ...)(S(I) + hC[e\ + ...)(/ + hr T + ...) + ... 

Since residual stress has been assumed away, this simplifies to 

= hC[e\ + ... 

which establishes the claim, i.e. that the linear part of the stress depends only on the 
symmetric part of Vw. 

This fact implies in turn that C^ki = (Problem 1). 

This symmetric part of Vu is called the infinitesimal strain: indeed, the Lagrangian 
strain is (Chapter 2) 

£ = i((Vxf(Vx)-J) 

= ^(/i(Vu + (Vuf) + /i 2 (Vuf(Vu)) 
= he + ... 

Define a corresponding infinitesimal stress by assuming an expansion of S in h: 

S = ha + ... 
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Then the linearized relation between stress and motion becomes 

a = C[e], i.e. = ^Cijki^ki 
ki 

which is often called Hooke's law: it simply states that, in the small motion limit, stress 
(force) is proportional to strain (relative displacement). 

To understand the dynamical laws (conservation of mass, momentum) in the small motion 
limit, it is necessary to compute the Taylor series of the velocity field in h: 

v(x,i) = — (q(x,t),f) 

h -fa(*- hu(q(x,t),t)) 

, <9u . 
= h—(x,t) + ... 

Note the simplicity of this result: to first order in h, the velocity is simply the time derivative 
of displacement! 

Introduce the expansion p(x, t) = p (x.,t) + hpi(x.,t] h) of the density. Then 

£f. + (V • v)p = ^ + v • Vp + (V • v)p 

Since this must hold identically in h, the coefficient of each power of h must vanish. The 
constant term (with respect to h) is 

£- 

and this is all that conservation of mass amounts to in the linear limit: as you will see 
shortly, only po plays a role in the conservation of momentum. So conservation of mass <^ 
mass density is time-independent to lowest order in h. 

Another way to see this is to recall that the spatial description of the conservation of 
mass is equivalent to the material description 

p(x, t) = Po(q(x, f)) det V x q(x, t) 

Expand x and its Jacobian in powers of h to see that I chose the right notation for p , i.e. 
p(x,t) = p (x) + .... 

As for conservation of momentum, 

=(» + •••) [ h W + -) 
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on , „ 

= hp ^ + ... = hV-a + ... + i 

at 2 

In order for this to work, clearly the external force f must vanish in the small motion limit, 
ie. f = M\ + .... Equating powers of h, obtain the linear elastodynamic equations of motion 

d 2 u 

Po ~df2 = V • o- + fi, a = C[e] 

It's customary to lose the subscripts on p and f , remembering that these quantities actually 
represent the lowest-order terms in h in the expansion of each. 

So far we have established these symmetries of the Hooke tensor: 

Cijki Cjiki 

Cijkl Cijlk 

These simply express the fact that C is a linear transformation on the space of symmetric 
3x3 matrices. Since the space of such matrices is 6 dimensional, the number of independent 
C's is 36 (the number of entries in a 6 x 6 matrix). [Idle amusement: convince yourself that 
this also follows from the two equations above.] 

Problems 5-10 show that conservation of energy, a concept introduced here for the first 
time in this course, implies additional symmetries: C is actually a symmetric linear trans- 
formation on Rgymm- 

Therefore there are actually 21 independent C^m (without any further assumptions about 
the nature of the elastic material). 



7.4 Isotropic Linear Elasticity 

Suppose that the elastic material undergoing small motion is isotropic. Expand both sides 
of the equation expressing isotropy in powers of h: 

RS(Vx)R T = hRC[Vu}R T + ... = hRC[e]R T + ... 

= S(RVxR T ) = hC[RVuR T ] + ... = hC[ReR T ] + ... 

Equating powers of h, and noting that e can be any symmetric matrix (why?), you see that 
C is an isotropic linear transformation on Rgy^ m : 

RC[e}R T = C[ReR T ] 
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for any rotation matrix R. We have already derived and used a structure theorem for 
isotropic transformations: C must take the form 

C[e\ = Atre/ + 2/ie 

for suitable scalars A and /i, called Lame parameters in this context. 

Problems 2 - 4 in this Chapter explain some of the physical significance of A and fi, 
from the point of view of statics - i.e. motions not dependent on t (often called simply 
"deformations"). The next sections give a complementary, dynamical point of view. 

7.5 Homogeneous isotropic media and alternate forms 
of the equations of motion 

Since 

tre = V • u 

and 

V • (//) = V/ 
for any scalar function /, you can rewrite 

V- (Atre/) = V(AV-u) 

If the material is homogeneous, i.e. neither p, A, nor p depend on x, then another simplifi- 
cation becomes possible: since 

V • Vu = V 2 u 

and 

V • (Vu) T = V(V • u) 
Thus for homogeneous isotropic elastic materials, the equations of motion become 

d 2 u 

P^ = (A + ^)V(V-u)+pV 2 u 

= (A + 2/i)V(V • u) - pV x V x u 
where in the last step we have used the ever-useful identity for vector valued functions w: 

V(V • w) = V 2 w + V x V x w 

(see Problem 4-6b). 
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7.6 Linear elastic waves in isotropic homogeneous ma- 
terials 



Motivated partly by the elastic fluid example, we hypothesize solutions of the plane wave 
form 

u(x,t) = f(ct-x-9)w 

where c is a positive scalar (the wave speed), # is a constant vector of unit length (the wave 
vector = normal to the equal phase surfaces), and w is a constant vector (the polarization 
vector = direction of particle motion). Note that 

V • (/w) = V/ • w = -f'6 • w 

V(V-(/w))=/"0(0.w) 
V x (/w) = -f'6 x w 

V x V x (/w) = f"9 x9xw 

and 

dt^ 1 

(where in each of these equations I have suppressed the argument ct — x • 9 of /). 

Therefore for functions of this form, the equations of motion are equivalent (after dividing 
through by p) to 

f(c 2 w) = /"( (^^j ee-w-^exexw) 

If /" does not vanish, then w must be an eigenvector with eigenvalue c 2 of the matrix of the 
linear operator appearing in the last equation. Note that 

w !->■ 9(9 ■ w) 

has the eigenvalue 1 with multiplicity 1 (vectors parallel to 0) and with multiplicity 2 
(vectors perpindicular to 9), whereas 

w h-> 9 x 9 x w) 

has eigenvalue -1 with multiplicity 2 (vectors perpindicular to 9) and with multiplicity 
1 (vectors parallel to 9). Since these matrices have the same eigenvectors, the eigenvalues 
of the linear combination of the matrices are simply the same linear combination of the 
eigenvalues. So we get two eigenvalues: 



• ^y^, multiplicity 1, corresponding to wave speed c = a = y— ^ and wave vector 9. 
These waves move in the same direction as the particle motion (i.e. the wave vector 
and the polarization vector are parallel), and are called compressional or longitudinal 
waves. 
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• ^, multiplicity 2, corresponding to wave speed c = /3 = and wave vectors per- 
pindicular to 9. These waves move in a direction perpindicular to particle motion (i.e 
the wave vector and the polarization vector are perpindicular) and are called shear or 
transverse waves. 



These solutions might seem to be quite special. However arbitrary solutions may be 
approximated arbitrarily well by linear combinations of these (Ch. 4), so they are in some 
sense not at all special! 



7.7 A few words about anisotropic and/or heteroge- 
neous materials 

Anisotropic (i.e. not isotropic) materials abound in nature, and in general require all 21 
possible independent Hooke's tensor coefficients to characterize them. Various partial sym- 
metries are possible: for example, shales (rocks composed of tiny platelets of mica, glued 
together with clay and other materials) have approximate rotational symmetry about the 
axis perpindicular to the bedding plane in which the platelets lie ("transverse isotropy"). 
Crystalline salt has the symmtries of a cubic lattice. These symmetries restrict the Hookc 
tensor in various complicated and interesting ways [refs]. 

Except where material homogeneity was explicitly assumed, nothing in the derivations 
in the section depends on it. So the density p (really po of course) and the Hooke tensor 
Cijki may depend on position x. 

The significance of the Lame parameters, being microscopic, is independent of homogene- 
ity also. However the construction of wave solutions in the last section depended explicitly 
on homogeneity, so the existence of elastic waves is not at this point established for hetero- 
geneous materials. 



7.8 Problems: 



Problem 7-1 .In the notes it was shown that, for an elastic material, the infinitesimal 
stress a depends only on the infinitesimal strain e. In terms of the Hooke's tensor C, this 
means that for any displacement u(x, t), 



dS 
dL 



(Vu) = C[Vu] = C 



1 



(Vu+(Vu) 7 ) 
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Show that (1) for any matrix L, there is a 
configuration so that V q x(q , 0) = I + L, 
displacement u(x, t), a time to, and a point 
L = Vu(x ,to)) and finally (3) therefore for 

C[L\ = C 

whence 

Cijki Cijik-, 



motion x(q, t) and a point q in the reference 
(2) consequently for any matrix L there is a 
x in the material volume at time to, for which 
any matrix L, 

> + ^) 

hJi k,l — 1, 3 



NOTE: In the next three problems, time plays no role, so I drop it from the notation. 



Problem 7-2 (shear modulus) . Define an displacement u by 

u(x) = (0,7^, Of, 7 G R 

Draw a picture of this deformation (see Problem 2-2 for a related picture!). Compute e. 
Assume that the material is an isotropic elastic solid, and show that in that case a = 2/xe. 

This deformation is infinitesimal shear. The Lame constant /i gives the relation between 
stress and strain for such deformations. Put another way: /i measures the force generated 
by, or equivalently the resistance to, a shear motion. So ji is called the shear modulus. 



Problem 7-3 (bulk modulus) . Define the extension displacement u of an isotropic 
elastic solid by 

u(x) = (7x1,7x2,7x3)^ 7 G R 
Show that e = 7/, a = 3k,€ where 

. 2 

K = \ + -fi 

Since k measures the resistance of the material to extension, which changes the volume 
without any shear deformation, it is called the bulk modulus of the material. 



Problem 7-4 (Poisson's ratio) . Show that that a pure tension infinitesimal stress in 
the x\ direction of an isotropic elastic solid, defined by 

a = diag(l,0,0) 

is the result of an infinitesimal deformation whose strain takes the form 

1 v v 



diag 



E E E 
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where 



E = 



/i(3A + 2/i) 
A + 11 



is Young's modulus, and 



A 



2(A + /i) 



is Poisson's ratio. That is, Young's modulus measures the increase in length resulting from 
a unit force along an axis, and Possion's ratio controls the decrease in diameter of the cross 
section perpindicular to the force axis. 

Young's modulus is always positive. Almost all materials have positive Poisson's ratio, 
i.e. they get skinnier if you pull on them. 

The next few problems lead to the relation between conservation of energy and symme- 
tries of the Hooke tensor Cijki for linearly elastic materials. Time is back in play. 

In linear elasticity, the potential energy of the material in Q C R 3 is 



(note that you can write the subscript i,j instead of j,i in the last factor because e is 
symmetric). Note that this quantity expresses the work done by the displacement (strain) 
against the internal forces (stress), and has the correct units for energy. 

The kinetic energy is (as usual) 



which should hold for a domain f2 in which there is no motion near the boundary, so that 
energy is not being transported out of the domain. That is, assume in the following that 
u = near dQ. 



Problem 7-5. Use calculus and Hooke's law to show that 





Conservation of (total mechanical) energy in this model is expressed as 



d_ 
~dt 



{KE + PE) = 
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— PE = - / tr — e + cr— 

] r . . r 



Problem 7-6. Use the symmetry of o and the definition of e to show that 

de / <9u \ 

trff S = tr<T v » 



Problem 7-7. Use the divergence theorem in some form to show that for any matrix 
valued function A, and any vector valued function w vanishing near the boundary dfl, 



( tr[(Vw)A] = - / w(V- A) 
Jn Jn 



Problem 7-8. Combine the results of the last several problems to conclude that 

Z PB = \J a nc», -/£■<?■*) 



ijkl 



Problem 7-9. Apply momentum balance to conclude that 

Therefore energy is conserved if the Hooke tensor Cijki has the additional symmetry 

Cijkl Cklij 
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Problem 7-10. Show that 

u(x,t) = (A + Bt)x 

is a motion consistent with momentum balance, ie. satisfying the equations of motion, for 
any symmetric matrices A and B. Show that conservation of energy then implies that 

y^ACkig — Cijki)A k iBij = 

ijkl 

and that since A and B are arbitrary symmetric matrices, this implies that 

Cijkl Cklij- 

Thus this symmetry is equivalent to conservation of energy. 
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Chapter 8 



Thermodynamics and the Heat 
Equation 



The goal of this chapter is to show how the heat equation arises from the first and second 
laws of thermodynamics and various contitutive assumptions about heat and deformation. 



as was discussed in Chapter 7. The potential energy however may comprise both strain 
energy as introduced there, even in a nonlinearized variant, and energy associated with the 
vibrational states of molecules. This latter motion is not represented by strain nor by any 
other component of continuum theory, so we must introduce a density e which is supposed 
to model it, together with any macroscopic (strain) contributions to energy. By convention, 
e has units of energy per unit mass, whence the total energy of u t is 



The power supplied to oo t is the rate of work done on it, and is supposed to consist of the 
work done on the material in u t by the material outside of it, and any work done by an 
otherwise unspecified body force represented by a force volume density f . A small piece of 
material centered at x moves approximately through the displacement Atv(x, t) between 
times t and t + At, and experiences the force f(x, t)AV, AV being the volume of the bit 



8.1 The First Law 



The kinetic energy of a moving part u t is 





85 



of material. Thus the work done on the material due to the motion and extenal force is 
approximately f(x, t) ■ v(x, t)AVAt. Adding up all such contributions and dividing by At, 
obtain the rate of work done on uj t by f: 

/f.v 

Jut 

Similarly, the force on a small area AA at x on the surface of Ut is AS'(x)n so by a similar 
argument the rate of work done on this small area is (Sn) ■ v, and on the entire moving part 
the surface integral of this quantity. 

The remaining potential energy tied up in the material is supposed to take the form 
of heat, for which there is both a volume source and a surface source. The surface source 
(i.e. heat transferred through the boundary of the moving part, due to contact with the 
remaining material) is supposed to be represented by an area density which is a function 
only of the surface normal. We supposed this for the traction (Chapter 3), and by a similar 
argument to that used by Cauchy conclude that the surface source of heat must take the 
form — h ■ n, in which h is the heat flux vectorfield. The sign convention means that h 
represents the direction of heat flow, whence its contribution to the energy in uj t is negative 
if the flux vector h points in the same direction as the outwoard normal. The volume source 
is supposed to represent processes tied to the material, for instance chemical reactions, is 
therefore a density r per unit mass. 

The upshot of these considerations is the integral form of the principle of conservation 
of energy, also called the first law of thermodynamics: the rate of change of energy in the 
moving part is the sum of the four rates described. In equations, 

4 / p(^IN 2 + e) = / (f-v + pr)+ / (Sn-h)-v 

A differential form follows from Reynold's Transport Theorem, esp. problem 1-6, and the 
Divergence Theorem in a manner that must by now seem familiar. The only novelty is an 
identity for the divergence of (Sn) ■ v = (Sv) ■ n (why?), which is 

V • (Sv) = (V • S) ■ v + S ■ Vv 

Here the second "•" is the Frobenius inner product of matrices 

A ■ B = tr A T B = AijBij 

We obtain 

/ Dv De\ 
P V' Dt + Dt) =(f + V-5)-v + 5-Vv + pr-V-h 

This simplifies still further, in view of the conservation of linear momentum, and we get the 
final differential form of the first law: 

Be 

p— = S ■ Vv + pr - V • h 
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8.2 The Second Law 



Thermodynamics posits the existence of a second fundamental material field, the entropy of 
the material, denoted 77. Very crudely speaking, entropy is meant to measure the relative 
disorder of the microscopic motion leading to the heat contribution to energy. It is sup- 
posed to have volume and surface sources in any moving part ou t , and furthermore these are 
proportional to the similar quantities for heat by the reciprocal of the temperature 9. The 
difference between the first and second laws is that the rate of entropy is not supposed to 
be equal to these source contributions, but at least equal to it: 



-/ pv> I --/ — 

dt JU)t J Lit 



/ojt Jojt v Jdojt 

For this to make sense, the temperature must be presumed positive. 

This inequality is the integral statement of the second law of thermodynamics, also known 
as the Clausius-Duhem inequality. 

After familiar manipulations, the second law takes the differential form: 

Dn V6 

1 Dt ~ pr + V ' h ~ h ' IT- 

Two of the terms on the right hand side also appear in the first law (conservation of energy), 
i.e. 

„ , _ De 

-pr + V • h = S ■ Vv - p— 

Combining these last two relations yields 

/ Dn De\ V6 

The term in parentheses on the left is (almost) the material derivative of the Helmholtz free 
energy ip = e — 9rj: 

Dili _ De Dr] DO 

which leads to another local form of the second law: 

Dip D6 n , W 



8.3 Stationary, undeformed heat conductor 

The relations introduced so far between the various material fields of the theory - v, S, e - are 
insufficient to predict their evolution, i.e. their values at t > t Q given their values at t — t . 
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A count of variables vs. equations suggests that additional relations are required, as was the 
case in all of the preceding chapters. Specification of these additional consitutive laws is a 
fairly complex task in general - for example, they must satisfy the principle of independence 
of observer, which strongly constrains their general form. The matter is a bit simpler in case 
the material undergoes no motion at all, and only the transfer of heat is to be perdicted. 

"No motion" means: v = 0, material and partial (time) derivatives coincide for all fields, 
material and spatial descriptions coincide, and all macroscopic mechanical contributions to 
energy balance vanish. Conservation of energy takes the form 

de _ , 

p- = pr-V-h 



and the second law may be written 



dip 09 , V0 n 



We make a major supposition: that any time "snapshot" of the fields 9, ip, and 77, for any 
fixed t may be reached at steady state, i.e. when these fields are in fact unchanging. This 
presumes that sources and sinks (r, h) can actually be prescribed somehow to keep these 
fields constant, and equal to any particular values they might attain in any other way. 

This is a big supposition, and it is not a priori obvious why it should be acceptable. It 
can be justified a posteriori for the system that will be introduced below. 

"Steady state" means 

89 dip drj 
dt = ~dt = di = 
This being so, the second law boils down to 

-h ■ V0 > 

(recalling that 9 is positive). 

The important point is that this inequality must be true regardless of whether the system 
is in steady state or not, because we have assumed that any temperature field (along with 
the other fields) can occur in steady state if it occurs at all. 

Together with the first law, this inequality implies that heat flows from higher temperature 
to lower temperature regions!. 

Indeed, suppose that u t = {x : 9(x,t) > 9 }, and that 9 is not stationary near the 
boundary of ui t . Then at every point on the boundary, n = — V0/||V0|| (note that the 
minus sign makes the normal outward, as is required by the divergence theorem!). Thus the 
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second law implies that h ■ n > 0, whence the term on the right hand side of the first law, 
expressing the transfer of heat through the boundary of oo t to the remainder of the material, 
is non-positive. Thus heat flows from hot places to cold places, consistent with everyday 
experience. 



8.4 Fourier's Law and the Heat Equation 



The preceding argument makes it reasonable to suppose that the heat flux h is a function 
of the temperature gradient, and possibly of the temperature and the position in the mate- 
rial: The simplest form which this dependence could take, consistent with independence of 
observer, was posited by Fourier: 

h = -k(0,x)V0 

Independence of observer requires that k = kl, in which k is the heat conductivity and 
k(9, x) > to ensure that — h • V9 > as required by the second law. 

Additional laws are required to determine tp and r\: the theory assumes that these are 
functions of 9 (and possibly of position). One could permit these quantities to depend on 
V9 as well, as does h, but independence of observer forces that dependence to be trivial. In 
view of the relation between ip, rj, and e, it follows that e = e(0,x), whence 

de _ 89 
dt~ C dt 

where c = 8e/89 is the specific heat of the material, and is assumed positive (rate of increase 
of heat energy with temperature). The specific heat has units of energy per unit mass per 
unit temperature. 

With these additional relations assumed, the first law reads 

39 

pc(9)— = V ■(k(9)V9)+pr 



For 9 varying by a small amount around an equilibrium value 6q, first order perturbation 
gives a simplified, linear PDE for 9. Assuming that (i) the reference temperature is in 
steady state (d9 /dt = 0) and (ii) the heat conductivity k(9) depends only weakly on 9 
(dk/89 negligible) you obtain 

89 

pc — = V-(k V9)+pr 

for the fluctuation 9 about 9 with c = c(9o),ko = k(9 ). Note that these quantities may 
still depend on x, which dependence accomodates variation in material properties from one 
part of the heat-conducting body to another. 
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This last equation is popularly known as the heat equation. Note that it takes exactly 
the same form as the dynamical law for density or pressure in single-phase porous flow. For 
a homogeneous body (p, c , k independent of x, it takes the even simpler form 

pc — = k \7 9 + pr 



8.5 Problems: 



Problem 9-1 . [30] Assume that cq — k — 1.0 and r = (homogeneous heat equation with 
normalized coefficients), and that space is d- dimensional (d — 1, 2, or 3 are the interesting 
cases). Show that 

solves 
for t > 0. 



Problem 9-2 . [40] Show that G(-,t) — > 5{x) as t — > + , i.e. for any continuous function 
4> : R d — >■ R vanishing outside of a large ball, 

lim^ ,t>o / ^ / G(x / ,t)0(x / ) = 0(0) 

More generally, show that 

lim t _> ,t>o / dx'G(x - x', t)0(x') = 0(x) 

[If you show the first identity, the second follows quickly by a simple change of variable!] 
You may use the following famous identity without proof, if you can see a use for it: 

/CO 2 
dxe~ x = \Jlx 
-co 



Problem 9-3 . [30] Show how the Gaussian kernel G should be modified so that the 
properties established in problems 9-1 and 9-2 continue to hold for general positive and 
constant (independent o/x, t) p,Co,k . 
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Problem 9-4 . [Extra Credit - 30] Show that the kernel constructed in problem 9-3 may- 
be used to solve the inhomogeneous heat equation as stated above: that is, that 

0(x, t) = /* dt' [ dx'G(x - x', t - t')pr{x.', t') + / dx'G(x - x', t)6 (x) 

JO JR" JR" 

solves the initial value problem for the heat equation 
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pc — = k W 2 e + pr, t > 0; 
at 

9(x,t) ->■ o (x), t^0,t>0 
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Chapter 9 



Porous Flow 



[Please see pp. 10-19 in Ewing's article "Problems arising in the modeling of processes for 
hydrocarbon recovery" which is Chapter 1 of The Mathematics of Reservoir Simulation, 
ed. R. Ewing, SIAM, 1983. I covered the material on single-phase flow in class; Ewing's 
discussion of multphase flow and chemical/compositional modeling is very readable. ] 
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Chapter 10 



Electrodynamics and Maxwell's 
Equations 



Maxwell's equations connect the following six fields: 

• electric, E(x, t) 

• magnetic, H(x, t) 

• electric, displacement D(x, t) 

• magnetic induction, B(x, t) 

• current density, j(x, t) 

• charge density, p(x, t) 

and involve a constant c with units of velocity ~ 3 x 10 5 km/s, which will be interpreted 
below as the speed of light in the vacuum. 

„ 19D 4tt. 

VxH--— = — j 

COT C 
C OT 

V • E = 4vrp; V • B = 

In addition several constitutive laws must be specified to fully determine the evolution of 
these fields. As in other topics in continuum mechanics, these express properties of particular 
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materials. For small intensity fields in isotropic materials moving slowly relative to c, they 
take the form 

j = o-E; D = eE; B = fiH 

in which the scalar fields a, e, fx are respectively the conductivity, the dielectric permittivity, 
and the magnetic permeability of the material. These satisfy the constraints a > 0, e, fx > 1. 
Electromagnetically speaking, a vacuum is a material in which these quantities take their 
extreme values: i.e. in a vacuum, a — 0, e = fx — 1 (in appropriate units). 

"See Jackson for details." 



10.1 Problems. 



Problem 1. [30] Suppose that all of the Maxwell fields are independent of time - as 
observed in class, this stationary assumption decouples the Maxwell equations governing 
the electric fields E, B from those governing the magnetic fields H, B. The two classes of 
fields remain linked however through the constitutive laws. Show that in this static situation, 
the electrostatic potential (f> defined by 

E = -V0 

satisfies 

V • (crV0) = 

and show that this relation can be interpreted as stating conservation of current in electro- 
statics: around any closed surface dfl bounding a domain Q, the net current circulating on 
the surface, defined as 

/ J ■ n, 

Jan 

vanishes. 



Problem 2. (i) Show that in a vacuum containing no charges and currents (j = 0,p = 0), 
the electric field satisfies the homogeneous wave equation 

1 * B - V>E = 0. 



c 2 dt 2 



Since we established in Ch. 5 that such equations have plane wave solutions moving with 
velocity c, the interpretation of c as the velocity of light (i.e. waves in the the electric field 
E) in the vacuum is at least partly justified. 
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(ii) Show how this equation should be modified in a homogeneous dielectric material, i.e. 
o" = 0. Note that the constraints on e, /i mentioned above imply that waves in such a material 
have velocity < c if the material is not a vacuum. 



Problem 3. [30] In class I defined a general scalar conservation law as a partial differential 
equation of the form 

du „ _ 
— + V-F = 

at 

and noted that we have encountered these throughout the semester, beginning with the law 
of conservation of mass. The divergence theorem implies that the rate of change of total u 
(whatever that may be) in a volume is the integral of F • n over the boundary, whence F is 
termed a flux. 

An interesting and much-studied special case is the nonlinear conservation law, in which 
F : R — > R 3 is a function of u. [Of course it could be a linear function, but that's not so 
interesting - and in fact is termed a "degenerate" case!]. 

(i) demonstrate the vector calculus identity 

OF 

V • (F(u(x, *))) = ^M x > *)) • V«(x, t) 



(ii) show that if the curve x : [0, T] — > R 3 is the solution of the ordinary differential equation 
then u is constant along this curve, i.e. w(x(t),t) = -u(x(0),0). 

(iii) deduce that any such curve (as defined in (ii)) is a straight line, and can be computed 
from the initial value u(x(0), 0). 

You can interpret this as a prescription for (more or less) explicitly solving the initial value 
problem for a nonlinear conservation law: given F and the initial data -u(x ,0), construct 
the straight line in spacetime 

(^ KXo ' 0)),t ) 

then w(x, t) = w(xo,0) at any point along this line. Such lines fill out spacetime, since you 
can start them anywhere on the t — surface. Therefore you have determined u. 

The fly in the ointment is that the lines may cross! Since this construction has to work 
if the solution u is smooth, the conclusion is that smooth solutions may not exist globally 
in time, even if the initial data is smooth. This observation leads to the concepts of weak 
solutions, shocks, and rarefaction waves, which will be discussed in CAAM 437. 
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